# -*- coding: utf-8 -*-
#
# Copyright (C) 2024 - 2025 ANSYS, Inc. and/or its affiliates.
# SPDX-License-Identifier: Apache-2.0
#
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from enum import Enum
import json
from pathlib import Path
import sys
import warnings
from ansys.aedt.core.aedt_logger import pyaedt_logger as logger
from ansys.aedt.core.generic.constants import AEDT_UNITS
from ansys.aedt.core.generic.constants import SpeedOfLight
from ansys.aedt.core.generic.constants import unit_converter
from ansys.aedt.core.generic.general_methods import conversion_function
from ansys.aedt.core.generic.general_methods import pyaedt_function_handler
from ansys.aedt.core.generic.numbers_utils import decompose_variable_value
from ansys.aedt.core.internal.checks import ERROR_GRAPHICS_REQUIRED
from ansys.aedt.core.internal.checks import check_graphics_available
from ansys.aedt.core.internal.checks import graphics_required
from ansys.aedt.core.visualization.plot.matplotlib import ReportPlotter
current_python_version = sys.version_info[:2]
if current_python_version < (3, 10): # pragma: no cover
raise Exception("Python 3.10 or later is required for radar postprocessing.")
try:
from scipy.interpolate import RegularGridInterpolator
except ImportError: # pragma: no cover
warnings.warn(
"The SciPy module is required to use the 'rcs_visualization.py' module.\nInstall with \n\npip install scipy"
)
try:
import numpy as np
except ImportError: # pragma: no cover
warnings.warn(
"The NumPy module is required to use the 'rcs_visualization.py' module.\nInstall with \n\npip install numpy"
)
np = None
# Check that graphics are available
try:
check_graphics_available()
import pyvista as pv
from ansys.tools.visualization_interface import MeshObjectPlot
from ansys.tools.visualization_interface import Plotter
from ansys.tools.visualization_interface.backends.pyvista import PyVistaBackend
except ImportError: # pragma: no cover
warnings.warn(ERROR_GRAPHICS_REQUIRED)
try:
import pandas as pd
except ImportError: # pragma: no cover
warnings.warn(
"The Pandas module is required to use the `rcs_visualization.py.` module. \nInstall with \n\npip install pandas"
)
pd = None
[docs]
class MonostaticRCSData(object):
"""Provides monostatic RCS data.
This class reads the monostatic RCS metadata in a JSON file
generated by the PyAEDT ``MonostaticRCSExporter`` class
and returns the Python interface to plot and analyze
the RCS data.
Parameters
----------
input_file : str
Metadata information in a JSON file.
Examples
--------
>>> from ansys.aedt.core import Hfss
>>> from ansys.aedt.toolkits.radar_explorer.rcs_visualization import MonostaticRCSData
>>> app = Hfss(version="2025.1", design="Antenna")
>>> data = app.get_rcs_data()
>>> metadata_file = data.metadata_file
>>> app.release_desktop()
>>> rcs_data = MonostaticRCSData(input_file=metadata_file)
"""
def __init__(self, input_file):
input_file = Path(input_file)
# Public
self.output_dir = input_file.parent
if not input_file.is_file():
raise FileNotFoundError("JSON file does not exist.")
# Private
self.__logger = logger
self.__input_file = input_file
self.__raw_data = {}
self.__frequency = None
self.__name = None
self.__solution = None
self.__incident_wave_theta = None
self.__incident_wave_phi = None
self.__available_incident_wave_theta = None
self.__available_incident_wave_phi = None
self.__frequencies = []
with input_file.open("r") as file:
self.__metadata = json.load(file)
self.__frequency_units = self.__metadata["frequency_units"]
self.__monostatic_file = None
if self.__metadata["monostatic_file"]:
self.__monostatic_file = self.output_dir / self.__metadata["monostatic_file"]
self.__data_conversion_function = "dB20"
self.__window = "Flat"
self.__window_size = 1024
self.__aspect_range = "Horizontal"
self.__upsample_range = 512
self.__upsample_azimuth = 64
self.__upsample_elevation = 64
self.__interpolation = "linear"
self.__extrapolate = True
self.__gridsize = "Middle"
if self.__monostatic_file and not self.__monostatic_file.is_file():
raise Exception("Monostatic file is invalid.")
self.rcs_column_names = ["data"]
# Load farfield data
if self.__monostatic_file:
is_rcs_loaded = self.__init_rcs()
else:
is_rcs_loaded = True
if not is_rcs_loaded: # pragma: no cover
raise RuntimeError("RCS information cannot be loaded.")
if self.__monostatic_file:
# Update active frequency if passed in the initialization
self.frequency = self.frequencies[0]
@property
def raw_data(self):
"""Antenna data."""
return self.__raw_data
@property
def metadata(self):
"""Antenna metadata."""
return self.__metadata
@property
def name(self):
"""Data name."""
if isinstance(self.raw_data, pd.DataFrame):
self.__name = self.raw_data.columns[0]
return self.__name
@property
def solution(self):
"""Data solution name."""
self.__solution = self.__metadata["solution"]
return self.__solution
@property
def input_file(self):
"""Input file."""
return self.__input_file
@property
def frequency_units(self):
"""Frequency units."""
return self.__frequency_units
@property
def frequencies(self):
"""Available frequencies."""
if isinstance(self.raw_data, pd.DataFrame) and "Freq" in self.raw_data.index.names:
frequencies = np.unique(np.array(self.raw_data.index.get_level_values("Freq")))
self.__frequencies = frequencies.tolist()
return self.__frequencies
@property
def available_incident_wave_theta(self):
"""Available incident wave theta."""
if isinstance(self.raw_data, pd.DataFrame) and "IWaveTheta" in self.raw_data.index.names:
self.__available_incident_wave_theta = np.unique(
np.array(self.raw_data.index.get_level_values("IWaveTheta"))
)
return self.__available_incident_wave_theta
@property
def incident_wave_theta(self):
"""Active incident wave theta."""
if isinstance(self.raw_data, pd.DataFrame) and self.__incident_wave_theta is None:
self.__incident_wave_theta = self.available_incident_wave_theta[0]
return self.__incident_wave_theta
@incident_wave_theta.setter
def incident_wave_theta(self, val):
"""Active incident wave theta."""
if val in self.available_incident_wave_theta:
self.__incident_wave_theta = val
else:
self.__logger.error("Value is not available.")
@property
def available_incident_wave_phi(self):
"""Available incident wave phi."""
if isinstance(self.raw_data, pd.DataFrame) and "IWavePhi" in self.raw_data.index.names:
self.__available_incident_wave_phi = np.unique(np.array(self.raw_data.index.get_level_values("IWavePhi")))
return self.__available_incident_wave_phi
@property
def incident_wave_phi(self):
"""Active incident wave phi."""
if isinstance(self.raw_data, pd.DataFrame) and self.__incident_wave_phi is None:
self.__incident_wave_phi = self.available_incident_wave_phi[0]
return self.__incident_wave_phi
@incident_wave_phi.setter
def incident_wave_phi(self, val):
"""Active incident wave phi."""
if val in self.available_incident_wave_phi:
self.__incident_wave_phi = val
else:
self.__logger.error("Value is not available.")
@property
def frequency(self):
"""Active frequency."""
return self.__frequency
@frequency.setter
def frequency(self, val):
if isinstance(val, str):
frequency, units = decompose_variable_value(val)
unit_converter(frequency, "Freq", units, self.frequency_units)
val = frequency
if val in self.frequencies:
self.__frequency = val
# self.__freq_index = self.frequencies.index(val)
else:
self.__logger.error("Frequency not available.")
@property
def data_conversion_function(self):
"""RCS data conversion function.
The available functions are:
- ``"dB10"``: Converts the data to decibels using the base 10 logarithm.
- ``"dB20"``: Converts the data to decibels using the base 20 logarithm.
- ``"abs"``: Computes the absolute value of the data.
- ``"real"``: Computes the real part of the data.
- ``"imag"``: Computes the imaginary part of the data.
- ``"norm"``: Normalizes the data to have values between 0 and 1.
- ``"ang"``: Computes the phase angle of the data in radians.
- ``"ang_deg"``: Computes the phase angle of the data in degrees.
"""
return self.__data_conversion_function
@data_conversion_function.setter
def data_conversion_function(self, val):
available_functions = ["dB10", "dB20", "abs", "real", "imag", "norm", "ang", "ang_deg", None]
if val in available_functions:
self.__data_conversion_function = val
@property
def interpolation(self):
"""Interpolation method."""
return self.__interpolation
@interpolation.setter
def interpolation(self, val):
if val in RegularGridInterpolator._ALL_METHODS:
self.__interpolation = val
else:
self.__logger.error("Interpolation method not available.")
@property
def extrapolate(self):
"""Extrapolation flag."""
return self.__extrapolate
@extrapolate.setter
def extrapolate(self, val):
if isinstance(val, bool):
self.__extrapolate = val
else:
self.__logger.error("Extrapolation flag must be a boolean value.")
@property
def gridsize(self):
"""Grid size for ISAR."""
return self.__gridsize
@gridsize.setter
def gridsize(self, val):
available_gridsizes = ["Inside", "Outside", "Middle"]
if val in available_gridsizes:
self.__gridsize = val
else:
self.__logger.error("Value for `gridsize` is invalid. The value must be 'Inside', 'Outside', or 'Middle'.")
@property
def window(self):
"""Window function.
The available functions are ``"Flat"``, ``"Hamming"``, and ``"Hann"``.
"""
return self.__window
@window.setter
def window(self, val):
available_functions = ["Flat", "Hamming", "Hann"]
if val in available_functions:
self.__window = val
else:
self.__logger.error("Value for `window` is invalid. The value must be 'Flat', 'Hamming', or 'Hann'.")
@property
def window_size(self):
"""Window size."""
return self.__window_size
@window_size.setter
def window_size(self, val):
self.__window_size = val
@property
def aspect_range(self):
"""Aspect range for ISAR."""
return self.__aspect_range
@aspect_range.setter
def aspect_range(self, val):
self.__aspect_range = val
@property
def upsample_range(self):
"""Upsample range for ISAR."""
return self.__upsample_range
@upsample_range.setter
def upsample_range(self, val):
self.__upsample_range = val
@property
def upsample_azimuth(self):
"""Upsample azimuth for ISAR."""
return self.__upsample_azimuth
@upsample_azimuth.setter
def upsample_azimuth(self, val):
self.__upsample_azimuth = val
@property
def upsample_elevation(self):
"""Upsample elevation for ISAR."""
return self.__upsample_elevation
@upsample_elevation.setter
def upsample_elevation(self, val):
self.__upsample_elevation = val
@property
def rcs(self):
"""RCS data for active frequency, theta, and phi."""
rcs_value = None
if isinstance(self.raw_data, pd.DataFrame):
data = self.rcs_active_theta_phi
data = data[data["Freq"] == self.frequency]
data = data.drop(columns=["Freq"])
value = data.values
rcs_value = value[0][0]
return rcs_value
@property
def rcs_active_theta_phi(self):
"""RCS data for active theta and phi."""
rcs_value = None
if isinstance(self.rcs_active_theta, pd.DataFrame):
data = self.rcs_active_theta
data = data[data["IWavePhi"] == self.incident_wave_phi]
rcs_value = data.drop(columns=["IWavePhi"])
return rcs_value
@property
def rcs_active_frequency(self):
"""RCS data for active frequency."""
value = None
if isinstance(self.raw_data, pd.DataFrame):
data = self.raw_data.xs(key=self.frequency, level="Freq")
data_converted = conversion_function(data[self.name], self.data_conversion_function)
df = data_converted.reset_index()
df.columns = ["IWavePhi", "IWaveTheta", "Data"]
value = df
return value
@property
def rcs_active_theta(self):
"""RCS data for active incident wave theta."""
value = None
if isinstance(self.raw_data, pd.DataFrame):
data = self.raw_data.xs(key=self.incident_wave_theta, level="IWaveTheta")
df = data.reset_index()
df.columns = ["Freq", "IWavePhi", "Data"]
value = df
return value
@property
def rcs_active_phi(self):
"""RCS data for active incident wave phi."""
value = None
if isinstance(self.raw_data, pd.DataFrame):
data = self.raw_data.xs(key=self.incident_wave_phi, level="IWavePhi")
df = data.reset_index()
df.columns = ["Freq", "IWaveTheta", "Data"]
value = df
return value
@property
def range_profile(self):
"""Range profile."""
value = None
if isinstance(self.raw_data, pd.DataFrame):
data = self.rcs_active_theta_phi["Data"]
# Take needed properties
size = self.window_size
nfreq = len(self.frequencies)
# Compute window
win_range, _ = self.window_function(self.window, nfreq)
windowed_data = data * win_range
# Perform FFT
sf_upsample = self.window_size / nfreq
windowed_data = np.fft.fftshift(sf_upsample * np.fft.ifft(windowed_data.to_numpy(), n=size))
data_converted = conversion_function(windowed_data, self.data_conversion_function)
df = unit_converter((self.frequencies[1] - self.frequencies[0]), "Freq", self.frequency_units, "Hz")
pd_t = 1.0 / df
dt = pd_t / size
range_norm = dt * np.linspace(start=-0.5 * size, stop=0.5 * size - 1, num=size) / 2 * SpeedOfLight
index_names = ["Range", "Data"]
df = pd.DataFrame(columns=index_names)
df["Range"] = range_norm
df["Data"] = data_converted
value = df
return value
@property
def waterfall(self):
"""Waterfall."""
waterfall_df = None
if isinstance(self.raw_data, pd.DataFrame):
waterfall_data = []
if self.aspect_range == "Horizontal":
original_phi = self.incident_wave_phi
for phi in self.available_incident_wave_phi:
self.incident_wave_phi = phi
range_profile = self.range_profile
new_data = {
"Range": range_profile["Range"],
"Data": range_profile["Data"],
"IWavePhi": np.full(range_profile["Range"].size, phi),
}
waterfall_data.append(pd.DataFrame(new_data))
self.incident_wave_phi = original_phi
else:
original_theta = self.incident_wave_theta
for theta in self.available_incident_wave_theta:
self.incident_wave_theta = theta
range_profile = self.range_profile
new_data = {
"Range": range_profile["Range"],
"Data": range_profile["Data"],
"IWaveTheta": np.full(range_profile["Range"].size, theta),
}
waterfall_data.append(pd.DataFrame(new_data))
self.incident_wave_theta = original_theta
waterfall_df = pd.concat(waterfall_data, ignore_index=True)
return waterfall_df
[docs]
@staticmethod
def phase_shift(data, slice_idx, num_after_fft, num_before_fft):
"""Apply phase shift to the data if the output has an even number of points in the ``slice_idx`` direction."""
if num_after_fft % 2 == 0:
d_arg = np.pi / num_after_fft
tmp = np.floor(-0.5 * num_before_fft) * d_arg
rng_shift_base = complex(np.cos(tmp), np.sin(tmp))
rng_shift_delta = complex(np.cos(d_arg), np.sin(d_arg))
for i in range(num_before_fft):
idx = [slice(None)] * data.ndim
idx[slice_idx] = i # this will pick [:, i] or [i, :] based on slice_idx
data[tuple(idx)] *= rng_shift_base
rng_shift_base *= rng_shift_delta
return data
@property
def isar_2d(self):
"""ISAR 2D."""
df = None
if isinstance(self.raw_data, pd.DataFrame):
# output size
ndrng = self.upsample_range
# get the input data
freqs = np.array(unit_converter(self.frequencies, "Freq", self.frequency_units, "Hz"))
nfreq = len(freqs)
freqs = freqs
avail_phi = np.array(self.available_incident_wave_phi)
avail_theta = np.array(self.available_incident_wave_theta)
# We are trying to create a 2D ISAR CUT potentially based off 3D ISAR data. When we choose
# a cut other than the main cuts, either for phi = 0 or theta = 90, we need to
# interpolate the data to the desired cut, which lies on a great arc. Given that potentially
# we'd have to extrapolate data over a large portion of the cut, we choose to disable
# extrapolation by default, and just zero pad, except if the data lies explicitly on one
# of the two main cuts. Because of numerical precision, even on the main cuts the resulting
# interpolants won't line up precisely with the original data, but we can forgive that for the sake of
# result accuracy and allow extrapolation.
if self.aspect_range == "Horizontal":
# Great arc (equatorial) cut, rotated by theta0 about y-axis
# We choose an arc with the same length and spacing as the original data, to preserve the
# original data's resolution and image extents.
el0 = np.deg2rad(90.0 - self.incident_wave_theta)
nangles = len(avail_phi)
arc_angles = np.deg2rad(avail_phi)
# Start with equatorial great circle (theta=90, phi=arc_angles)
x0 = np.cos(arc_angles)
y0 = np.sin(arc_angles)
# Rotate by el0 about y-axis
x = x0 * np.cos(el0)
y = y0
z = x0 * np.sin(el0)
# Convert to (theta, phi)
thetas = np.arccos(np.clip(z, -1.0, 1.0))
phis = np.arctan2(y, x)
theta_deg = np.rad2deg(thetas)
phi_deg = np.rad2deg(phis)
interp = RegularGridInterpolator(
(freqs, avail_phi, avail_theta),
self.raw_data.to_numpy().reshape((nfreq, len(avail_phi), len(avail_theta))),
method=self.interpolation,
bounds_error=False,
fill_value=None if self.incident_wave_theta == 90.0 else 0.0,
)
azel_samples = -arc_angles # convert from phi to azimuth
nxrng = nangles
else:
# Great arc (equatorial) cut logic for elevation
# We choose an arc with the same length and spacing as the original data, to preserve the
# original data's resolution and image extents.
phi0 = np.deg2rad(self.incident_wave_phi)
nangles = len(avail_theta)
arc_angles = np.deg2rad(avail_theta)
# Start with equatorial great circle (phi=0, theta=arc_angles)
x0 = np.sin(arc_angles)
z0 = np.cos(arc_angles)
# Rotate by phi0 about z-axis
x = x0 * np.cos(phi0)
y = x0 * np.sin(phi0)
z = z0
# Convert to (theta, phi)
thetas = np.arccos(np.clip(z, -1.0, 1.0))
phis = np.arctan2(y, x)
theta_deg = np.rad2deg(thetas)
phi_deg = np.rad2deg(phis)
interp = RegularGridInterpolator(
(freqs, avail_phi, avail_theta),
self.raw_data.to_numpy().reshape((nfreq, len(avail_phi), len(avail_theta))),
method=self.interpolation,
bounds_error=False,
fill_value=None if self.incident_wave_phi == 0.0 else 0.0,
)
azel_samples = np.pi / 2 - arc_angles # convert from theta to elevation
nxrng = nangles
# Interpolate for all frequencies over the points in the line el_deg, az_deg to obtain
# data along the desired cut.
data = np.empty((nfreq, nangles), dtype=np.complex128)
for i, freq in enumerate(freqs):
freq_points = np.column_stack([np.full_like(theta_deg, freq), phi_deg, theta_deg])
data[i, :] = interp(freq_points)
# center the cut samples, just in case
azel_samples = np.unwrap(azel_samples)
azel_samples -= np.mean(azel_samples)
fxmin = np.min(freqs)
fxmax = np.max(freqs)
if self.gridsize == "Inside":
f_ref = freqs[0]
elif self.gridsize == "Outside":
f_ref = freqs[-1]
else: # self.gridsize == "Middle" which matches the extents
f_ref = freqs[len(freqs) // 2]
fymax = np.sin(np.max(azel_samples)) * f_ref
fymin = np.sin(np.min(azel_samples)) * f_ref
fx = np.linspace(fxmin, fxmax, nfreq) # desired downrange frequencies
fy = np.linspace(fymin, fymax, nangles) # desired crossrange frequencies
f_x, f_y = np.meshgrid(fx, fy, indexing="ij")
# convert to equivalent freq and azimuthinterpolation points,
# so we can use the regular interpolator directly from freqs, az
fi = np.sqrt(f_x**2 + f_y**2)
azi = np.arctan2(f_y, f_x)
interp = RegularGridInterpolator(
(freqs, azel_samples),
data,
method=self.interpolation,
bounds_error=False,
fill_value=None if self.extrapolate else 0.0,
)
rdata = interp((fi, azi))
# Zero padding
if ndrng < nfreq:
# Warning('nx should be at least as large as the length of f -- increasing nx')
self.__logger.warning("nx should be at least as large as the number of frequencies.")
ndrng = nfreq
if nxrng < nangles: # pragma: no cover
# warning('ny should be at least as large as the length of az -- increasing ny');
self.__logger.warning("ny should be at least as large as the number of azimuth angles.")
nxrng = nangles
# Compute the image plane downrange and cross-range distance vectors (in
# meters)
dfx = fx[1] - fx[0] # difference in x-frequencies
dfy = fy[1] - fy[0] # difference in y-frequencies
dx = SpeedOfLight / (2 * dfx) / ndrng
dy = SpeedOfLight / (2 * dfy) / nxrng
x = np.linspace(start=0, stop=ndrng * dx, num=ndrng) # ndrng
y = np.linspace(start=0, stop=nxrng * dy, num=nxrng) # nxrng
# We want the physical extents of the image to be centered at the global origin, because
# that's how we draw the extents of the 2D ISAR domain.
# The center of the first pixel in the second half of each domain is centered at zero
# if the domain has odd length, but it is at dx/2 otherwise.
fft_domain_sizes = ((ndrng, nfreq), (nxrng, nangles))
for islide, fft_domain_size in enumerate(fft_domain_sizes):
rdata = MonostaticRCSData.phase_shift(rdata, islide, *fft_domain_size)
winx, winx_sum = self.window_function(self.window, nfreq)
winy, winy_sum = self.window_function(self.window, nangles)
winx = winx.reshape(-1, 1)
winy = winy.reshape(1, -1)
iq = np.zeros((ndrng, nxrng), dtype=np.complex128)
xshift = (ndrng - nfreq) // 2
yshift = (nxrng - nangles) // 2
iq[xshift : xshift + nfreq, yshift : yshift + nangles] = np.multiply(rdata, winx * winy)
# Normalize so that unit amplitude scatterers have about unit amplitude in
# the image (normalized for the windows). The "about" comes because of the
# truncation of the polar shape into a rectangular shape.
#
iq = np.fft.fftshift(iq) * ndrng * nxrng / winx_sum / winy_sum
isar_image = np.fft.fftshift(np.fft.ifft2(iq))
# Nx x Ny
isar_image = conversion_function(isar_image, self.data_conversion_function)
# isar_image = isar_image.transpose()
# isar_image = isar_image[::-1, :] # this used to be flipped, but it matched range/cross range defs now
# bring the center of the PHYSICAL image to 0, which means the first pixel on the
# second half is not at 0 for even length domains
range_values = x - 0.5 * (x[-1] - x[0])
cross_range_values = y - 0.5 * (y[-1] - y[0])
rr, xr = np.meshgrid(range_values, cross_range_values)
rr_flat = rr.ravel()
xr_flat = xr.ravel()
isar_image_flat = isar_image.ravel()
index_names = ["Down-range", "Cross-range", "Data"]
df = pd.DataFrame(columns=index_names)
df["Down-range"] = rr_flat
df["Cross-range"] = xr_flat
df["Data"] = isar_image_flat
return df
@property
def isar_3d(self):
"""ISAR 3D."""
df = None
if isinstance(self.raw_data, pd.DataFrame):
# get the input data
freqs = np.array(unit_converter(self.frequencies, "Freq", self.frequency_units, "Hz"))
nfreq = len(freqs)
az = np.unwrap(np.radians(-self.available_incident_wave_phi))
az -= np.mean(az)
nphi = len(az)
el = np.unwrap(np.radians(90.0 - self.available_incident_wave_theta))
el -= np.mean(el)
ntheta = len(el)
data = self.__raw_data
# output size
ndrng = self.upsample_range
nxrng1 = self.upsample_azimuth
nxrng2 = self.upsample_elevation
fxmin = np.min(freqs)
fxmax = np.max(freqs)
if self.gridsize == "Inside":
f_ref = freqs[0]
elif self.gridsize == "Outside":
f_ref = freqs[-1]
else: # self.gridsize == "Middle" which matches the extents
f_ref = freqs[len(freqs) // 2]
fymin = np.sin(np.min(az)) * f_ref
fymax = np.sin(np.max(az)) * f_ref
fzmin = np.sin(np.min(el)) * f_ref
fzmax = np.sin(np.max(el)) * f_ref
fx = np.linspace(fxmin, fxmax, nfreq) # desired downrange frequencies
fy = np.linspace(fymin, fymax, nphi) # desired crossrange frequencies
fz = np.linspace(fzmin, fzmax, ntheta) # desired crossrange frequencies
f_x, f_y, f_z = np.meshgrid(fx, fy, fz, indexing="ij")
# convert to equivalent freq, azimuth, and elevation interpolation points,
# so we can use the regular interpolator directly from freqs, az and el
fi = np.sqrt(f_x**2 + f_y**2 + f_z**2)
azi = np.arctan2(f_y, f_x)
with np.errstate(invalid="ignore"):
eli = np.arcsin(np.where(fi != 0, -f_z / fi, 0.0))
interp = RegularGridInterpolator(
(freqs, az, el),
data.to_numpy().reshape((nfreq, nphi, ntheta)),
method=self.interpolation,
bounds_error=False,
fill_value=None if self.extrapolate else 0.0,
)
rdata = interp((fi, azi, eli))
if ndrng < nfreq:
# Warning('nx should be at least as large as the length of f -- increasing nx')
self.__logger.warning("nx should be at least as large as the number of frequencies.")
ndrng = nfreq
if nxrng1 < nphi:
# warning('ny should be at least as large as the length of az -- increasing ny');
self.__logger.warning("ny should be at least as large as the number of azimuth angles.")
nxrng1 = nphi
if nxrng2 < ntheta:
# warning('nz should be at least as large as the length of el -- increasing nz');
self.__logger.warning("nz should be at least as large as the number of elevation angles.")
nxrng2 = ntheta
# We want the physical extents of the image to be centered at the global origin, because
# that's how we draw the extents of the 3D ISAR domain.
# The center of the first pixel in the second half of each domain is centered at zero
# if the domain has odd length, but it is at dx/2 otherwise.
fft_domain_sizes = ((ndrng, nfreq), (nxrng1, nphi), (nxrng2, ntheta))
for islide, fft_domain_size in enumerate(fft_domain_sizes):
rdata = MonostaticRCSData.phase_shift(rdata, islide, *fft_domain_size)
#
# add windowing
#
winx, winx_sum = self.window_function(self.window, nfreq)
winy, winy_sum = self.window_function(self.window, nphi)
winz, winz_sum = self.window_function(self.window, ntheta)
winx = winx.reshape(-1, 1, 1)
winy = winy.reshape(1, -1, 1)
winz = winz.reshape(1, 1, -1)
iq = np.zeros((ndrng, nxrng1, nxrng2), dtype=np.complex128)
xshift = (ndrng - nfreq) // 2
yshift = (nxrng1 - nphi) // 2
zshift = (nxrng2 - ntheta) // 2
window = winx * winy * winz
iq[xshift : xshift + nfreq, yshift : yshift + nphi, zshift : zshift + ntheta] = np.multiply(rdata, window)
#
# normalize so that unit amplitude scatterers have about unit amplitude in
# the image (normalized for the windows). The "about" comes because of the
# truncation of the polar shape into a rectangular shape.
#
iq = np.fft.fftshift(iq) * ndrng * nxrng1 * nxrng2 / winx_sum / winy_sum / winz_sum
isar_image = np.fft.fftshift(np.fft.ifftn(iq)) # Nx x Ny x Nz
isar_image = conversion_function(isar_image, self.data_conversion_function)
#
# compute the image plane downrange and crossrange distance vectors (in
# meters)
#
dfx = fx[1] - fx[0] # difference in x-frequencies
dfy = fy[1] - fy[0] # difference in y-frequencies
dfz = fz[1] - fz[0] # difference in z-frequencies
dx = SpeedOfLight / (2 * dfx) / ndrng
dy = SpeedOfLight / (2 * dfy) / nxrng1
dz = SpeedOfLight / (2 * dfz) / nxrng2
x = np.linspace(start=0, stop=ndrng * dx, num=ndrng)
y = np.linspace(start=0, stop=nxrng1 * dy, num=nxrng1)
z = np.linspace(start=0, stop=nxrng2 * dz, num=nxrng2)
#
# bring the center of the PHYSICAL image to 0, which means the first pixel on the
# second half is not at 0 for even length domains
range_values = x - 0.5 * (x[-1] - x[0])
cross_range1_values = y - 0.5 * (y[-1] - y[0])
cross_range2_values = z - 0.5 * (z[-1] - z[0])
rr, xr1, xr2 = np.meshgrid(range_values, cross_range1_values, cross_range2_values)
rr_flat = rr.ravel()
xr1_flat = xr1.ravel()
xr2_flat = xr2.ravel()
isar_image_flat = isar_image.ravel()
index_names = ["Down-range", "Cross-range-az", "Cross-range-el", "Data"]
df = pd.DataFrame(columns=index_names)
df["Down-range"] = rr_flat
df["Cross-range-az"] = xr1_flat
df["Cross-range-el"] = xr2_flat
df["Data"] = isar_image_flat
return df
[docs]
@staticmethod
def window_function(window="Flat", size=512):
"""Apply a window function.
Parameters
----------
window : str, default: ``"Flat"``
Unnormalized window function. Options are ``"Flat"``, ``"Hamming"``, and ``"Hann"``.
size : int, default: ``512``
Window size.
Returns
-------
tuple
Data windowed and data sum.
"""
if window == "Hann":
win = np.hanning(size)
elif window == "Hamming":
win = np.hamming(size)
else:
win = np.ones(size)
win_sum = np.sum(win)
return win, win_sum
@pyaedt_function_handler()
def __init_rcs(self):
"""Load monostatic radar cross-section data.
Returns
-------
bool
``True`` when successful, ``False`` when failed.
"""
try:
self.__raw_data = pd.read_hdf(self.__monostatic_file, key="df", mode="r")
except ImportError as e: # pragma: no cover
self.__logger.error(f"Failed to load monostatic RCS data: {e}")
return False
return True
[docs]
class MonostaticRCSPlotter(object):
"""Provides monostatic RCS plot functionalities.
Parameters
----------
rcs_data : :class:`ansys.aedt.toolkits.radar_explorer.rcs_visualization`, default: ``None``
Monostatic RCS data object.
Examples
--------
>>> from ansys.aedt.core import Hfss
>>> from ansys.aedt.toolkits.radar_explorer.rcs_visualization import MonostaticRCSData
>>> from ansys.aedt.toolkits.radar_explorer.rcs_visualization import MonostaticRCSPlotter
>>> app = Hfss(version="2025.1", design="Antenna")
>>> data = app.get_rcs_data()
>>> metadata_file = data.metadata_file
>>> app.release_desktop()
>>> rcs_data = MonostaticRCSData(input_file=metadata_file)
>>> rcs_plotter = MonostaticRCSPlotter(rcs_data)
"""
def __init__(self, rcs_data: MonostaticRCSData | None = None):
# Private
self.__rcs_data = rcs_data
self.__logger = logger
self.__model_units = "meter"
# Scene properties
self.show_geometry = True
self.__all_scene_actors = {"model": {}, "annotations": {}, "results": {}}
self.__x_max, self.__x_min, self.__y_max, self.__y_min, self.__z_max, self.__z_min = 1, -1, 1, -1, 1, -1
self.__model_info = None
self.__num_contours = 10
# Get geometries
if self.__rcs_data and "model_info" in self.rcs_data.metadata.keys():
self.__model_info = self.rcs_data.metadata["model_info"]
obj_meshes = self.__get_geometry()
self.__all_scene_actors["model"] = obj_meshes
# Get model extent
self.__get_model_extent()
@property
def rcs_data(self):
"""RCS data object."""
return self.__rcs_data
@property
def model_info(self):
"""Geometry information."""
return self.__model_info
@property
def model_units(self):
"""Model units."""
return self.__model_units
@property
def all_scene_actors(self):
"""All scene actors."""
return self.__all_scene_actors
@property
def extents(self):
"""Geometry extents."""
return [self.__x_min, self.__x_max, self.__y_min, self.__y_max, self.__z_min, self.__z_max]
@property
def center(self):
"""Geometry center."""
x_mid = (self.__x_max + self.__x_min) / 2
z_mid = (self.__z_max + self.__z_min) / 2
y_mid = (self.__y_max + self.__y_min) / 2
return np.array([x_mid, y_mid, z_mid])
@property
def radius(self):
"""Geometry radius."""
return max(
[abs(a) for a in (self.__x_min, self.__x_max, self.__y_min, self.__y_max, self.__z_min, self.__z_max)]
)
@property
def num_contours(self) -> int:
"""Number of contours."""
return self.__num_contours
@num_contours.setter
def num_contours(self, value):
self.__num_contours = value
[docs]
@pyaedt_function_handler()
def plot_rcs(
self,
primary_sweep="IWavePhi",
secondary_sweep="IWaveTheta",
secondary_sweep_value=None,
title="Monostatic RCS",
output_file=None,
show=True,
is_polar=False,
show_legend=True,
size=(1920, 1440),
):
"""Create a 2D plot of the monostatic RCS.
Parameters
----------
primary_sweep : str, default: ``"IWavePhi"``
X-axis variable. Options are ``"Freq"``, ``"IWavePhi"`` and ``"IWaveTheta"``.
secondary_sweep : str, default: ``"IWavePhi"``
Y-axis variable. Options are ``"Freq"``, ``"IWavePhi"`` and ``"IWaveTheta"``.
secondary_sweep_value : float, list, string, default: ``0``
List of cuts on the secondary sweep to plot. Options are
`"all"`, a single value float, or a list of float values.
title : str, default: ``"RectangularPlot"``
Plot title.
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image in not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
is_polar : bool, default: ``True``
Whether this plot is a polar plot.
show_legend : bool, default: ``True``
Whether to display the legend.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
Returns
-------
:class:`ansys.aedt.core.visualization.plot.matplotlib.ReportPlotter`
PyAEDT Matplotlib figure object.
"""
curves = []
all_secondary_sweep_value = secondary_sweep_value
if primary_sweep.casefold() == "freq" or primary_sweep.casefold() == "frequency":
x_key = "Freq"
x = self.rcs_data.frequencies
if secondary_sweep == "IWaveTheta":
if all_secondary_sweep_value is None:
all_secondary_sweep_value = self.rcs_data.incident_wave_theta
data = self.rcs_data.rcs_active_phi
data["Data"] = conversion_function(data["Data"], self.rcs_data.data_conversion_function)
y_key = "IWaveTheta"
else:
if all_secondary_sweep_value is None:
all_secondary_sweep_value = self.rcs_data.incident_wave_phi
data = self.rcs_data.rcs_active_theta
data["Data"] = conversion_function(data["Data"], self.rcs_data.data_conversion_function)
y_key = "IWavePhi"
else:
data = self.rcs_data.rcs_active_frequency
if primary_sweep.casefold() == "iwavephi":
x_key = "IWavePhi"
y_key = "IWaveTheta"
x = self.rcs_data.available_incident_wave_phi
if isinstance(secondary_sweep_value, str) and secondary_sweep_value == "all":
all_secondary_sweep_value = self.rcs_data.available_incident_wave_theta
elif secondary_sweep_value is None:
all_secondary_sweep_value = self.rcs_data.incident_wave_theta
else:
x_key = "IWaveTheta"
y_key = "IWavePhi"
x = self.rcs_data.available_incident_wave_theta
if isinstance(secondary_sweep_value, str) and secondary_sweep_value == "all":
all_secondary_sweep_value = self.rcs_data.available_incident_wave_phi
elif secondary_sweep_value is None:
all_secondary_sweep_value = self.rcs_data.incident_wave_phi
if all_secondary_sweep_value is not None:
if not isinstance(all_secondary_sweep_value, np.ndarray) and not isinstance(
all_secondary_sweep_value, list
):
all_secondary_sweep_value = [all_secondary_sweep_value]
for el in all_secondary_sweep_value:
data_sweep = data[data[y_key] == el]["Data"]
y = data_sweep.values
if not isinstance(y, np.ndarray): # pragma: no cover
raise Exception("Format of quantity is wrong.")
curves.append([x, y, "{}={}".format(y_key, el)])
if curves is not None:
new = ReportPlotter()
new.show_legend = show_legend
new.title = title
new.size = size
if is_polar:
props = {"x_label": x_key, "y_label": "RCS \n"}
for pdata in curves:
name = pdata[2] if len(pdata) > 2 else "Trace"
new.add_trace(pdata[:2], 0, props, name=name)
_ = new.plot_polar(traces=None, snapshot_path=output_file, show=show)
else:
from ansys.aedt.core.generic.constants import CSS4_COLORS
k = 0
for data in curves:
props = {"x_label": x_key, "y_label": "RCS", "line_color": list(CSS4_COLORS.keys())[k]}
k += 1
if k == len(list(CSS4_COLORS.keys())): # pragma: no cover
k = 0
name = data[2] if len(data) > 2 else "Trace"
new.add_trace(data[:2], 0, props, name)
_ = new.plot_2d(None, output_file, show)
return new
[docs]
@pyaedt_function_handler()
def plot_rcs_3d(self, title="Monostatic RCS 3D", output_file=None, show=True, size=(1920, 1440)):
"""Create a 3D plot of the monostatic RCS.
Parameters
----------
title : str, default: ``"RectangularPlot"``
Plot title.
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image is not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
Returns
-------
:class:`matplotlib.pyplot.Figure`
Matplotlib figure object.
If ``show=True``, a Matplotlib figure instance of the plot is returned.
If ``show=False``, the plotted curve is returned.
"""
data = self.rcs_data.rcs_active_frequency
rcs = data["Data"]
rcs_max = np.max(rcs)
rcs_min = np.min(rcs)
rcs_renorm = rcs + np.abs(rcs_min) if rcs_min else rcs
rcs_renorm = rcs_renorm.to_numpy()
# Negative values are not valid, this is cleaning numerical issues
rcs_renorm[rcs_renorm < 0] = 0.0
theta = np.deg2rad(data["IWaveTheta"])
phi = np.deg2rad(data["IWavePhi"])
unique_phi, unique_theta = np.unique(phi), np.unique(theta)
phi_grid, theta_grid = np.meshgrid(unique_phi, unique_theta)
r = np.full(phi_grid.shape, np.nan)
idx_theta = np.digitize(theta, unique_theta) - 1
idx_phi = np.digitize(phi, unique_phi) - 1
r[idx_theta, idx_phi] = rcs_renorm
x = r * np.sin(theta_grid) * np.cos(phi_grid)
y = r * np.sin(theta_grid) * np.sin(phi_grid)
z = r * np.cos(theta_grid)
new = ReportPlotter()
new.show_legend = True
new.title = title
new.size = size
quantity = f"RCS {self.rcs_data.data_conversion_function}"
props = {"x_label": "IWaveTheta", "y_label": "IWavePhi", "z_label": quantity}
new.add_trace([x, y, z], 2, props, title)
_ = new.plot_3d(trace=0, snapshot_path=output_file, show=show, color_map_limits=[rcs_min, rcs_max])
return new
[docs]
@pyaedt_function_handler()
def plot_range_profile(
self,
title="Range profile",
output_file=None,
show=True,
show_legend=True,
size=(1920, 1440),
):
"""Create a 2D plot of the range profile.
Parameters
----------
title : str, default: ``"RectangularPlot"``
Plot title.
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image is not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
show_legend : bool, default: ``True``
Whether to display the legend.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
Returns
-------
:class:`ansys.aedt.core.visualization.plot.matplotlib.ReportPlotter`
PyAEDT matplotlib figure object.
"""
data_range_profile = self.rcs_data.range_profile
ranges = np.unique(data_range_profile["Range"])
phi = self.rcs_data.incident_wave_phi
theta = self.rcs_data.incident_wave_theta
y = data_range_profile["Data"].to_numpy()
legend = f"Phi={np.round(phi, 3)} Theta={np.round(theta, 3)}"
curve = [ranges.tolist(), y.tolist(), legend]
new = ReportPlotter()
new.show_legend = show_legend
new.title = title
new.size = size
props = {"x_label": "Range (m)", "y_label": f"Range Profile ({self.rcs_data.data_conversion_function})"}
name = curve[2]
new.add_trace(curve[:2], 0, props, name)
_ = new.plot_2d(None, output_file, show)
return new
[docs]
@pyaedt_function_handler()
def plot_waterfall(
self, title="Waterfall", output_file=None, show=True, is_polar=False, size=(1920, 1440), figure=None
):
"""Create a 2D contour plot of the waterfall.
Parameters
----------
title : str, default: ``"RectangularPlot"``
Plot title.
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image is not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
is_polar : bool, default: ``True``
Whether to display in polar coordinates.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
figure : :class:`matplotlib.pyplot.Figure`, default: ``None``
An existing Matplotlib figure to add the plot to.
If no figure is provided, new ``Figure`` and ``Axes`` objects are created.
Returns
-------
:class:`ansys.aedt.core.visualization.plot.matplotlib.ReportPlotter`
PyAEDT Matplotlib figure object.
"""
data_range_waterfall = self.rcs_data.waterfall
ranges = np.unique(data_range_waterfall["Range"])
if self.rcs_data.aspect_range == "Horizontal":
phis = np.unique(data_range_waterfall["IWavePhi"])
ylabel = "Phi (deg)"
else:
phis = np.unique(data_range_waterfall["IWaveTheta"])
ylabel = "Theta (deg)"
phis = np.deg2rad(phis.tolist()) if is_polar else phis
n_range = len(ranges)
n_phi = len(phis)
values = data_range_waterfall["Data"].to_numpy()
values = values.reshape((n_range, n_phi), order="F")
ra, ph = np.meshgrid(ranges, phis)
if is_polar:
x = ph
y = ra
xlabel = " "
ylabel = " "
else:
x = ra
y = ph
xlabel = "Range (m)"
plot_data = [values.T, y, x]
new = ReportPlotter()
new.size = size
new.show_legend = False
new.title = title
props = {
"x_label": xlabel,
"y_label": ylabel,
}
new.add_trace(plot_data, 0, props)
_ = new.plot_contour(
trace=0,
polar=is_polar,
snapshot_path=output_file,
show=show,
figure=figure,
is_spherical=False,
max_theta=ra.max(),
min_theta=ra.min(),
)
return new
[docs]
@pyaedt_function_handler()
def plot_isar_2d(self, title="ISAR", output_file=None, show=True, size=(1920, 1440), figure=None):
"""Create a 2D pcolor plot of the 2D ISAR image.
Parameters
----------
title : str, default: ``"ISAR"``
Plot title.
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image is not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
figure : :class:`matplotlib.pyplot.Figure`, default: ``None``
An existing Matplotlib figure to add the plot to.
If no figure is provided, new ``Figure`` and ``Axes`` objects are created.
Returns
-------
:class:`ansys.aedt.core.visualization.plot.matplotlib.ReportPlotter`
PyAEDT matplotlib figure object.
"""
data_isar = self.rcs_data.isar_2d
ranges = np.unique(data_isar["Down-range"])
phis = np.unique(data_isar["Cross-range"])
n_range = len(ranges)
n_phi = len(phis)
values = data_isar["Data"].to_numpy().reshape((n_range, n_phi))
x, y = np.meshgrid(phis, ranges) # important, do not use ij here
xlabel = "Range (m)"
ylabel = "Cross Range (m)"
plot_data = [values, x, y]
new = ReportPlotter()
new.size = size
new.show_legend = False
new.title = title
props = {
"x_label": xlabel,
"y_label": ylabel,
}
new.add_trace(plot_data, 0, props)
_ = new.plot_pcolor(trace=0, snapshot_path=output_file, show=show, figure=figure)
return new
[docs]
@pyaedt_function_handler()
def plot_isar_3d(
self,
plane_cut: str = "xz",
plane_offset: float = 0,
title="ISAR 3D cut",
output_file=None,
show=True,
size=(1920, 1440),
figure=None,
):
"""Create a 2D pcolor plot of a 3D-ISAR cut.
Parameters
----------
title : str, default: ``"ISAR 3D cut"``
Plot title.
plane_cut : str, default: ``"xz"``
Cut to be plotted. Options are ``"xz"``, ``"yz"``, and ``"xy"``.
plane_offset : float, default: ``0``
Offset of the plane to plot. (It finds the closest available.)
output_file : str, default: ``None``
Full path for the image file. The default is ``None``, in which case an image in not exported.
show : bool, default: ``True``
Whether to show the plot. If ``False``, the Matplotlib instance of the plot is shown.
size : tuple, default: (1920, 1440)
Image size in pixel (width, height).
figure : :class:`matplotlib.pyplot.Figure`, default: ``None``
An existing Matplotlib figure to add the plot to.
If no figure is provided, new ``Figure`` and ``Axes`` objects are created.
Returns
-------
:class:`ansys.aedt.core.visualization.plot.matplotlib.ReportPlotter`
PyAEDT matplotlib figure object.
"""
# implementation for plane cut plot
if plane_cut is None:
raise ValueError("Please provide a plane cut type.")
if plane_offset is None:
raise ValueError("Please provide a plane offset value.")
data_isar_3d = self.rcs_data.isar_3d
down_range = data_isar_3d["Down-range"].unique()
cross_range_az = data_isar_3d["Cross-range-az"].unique()
cross_range_el = data_isar_3d["Cross-range-el"].unique()
# The above pivot may need to be adjusted depending on the data layout
# For now, let's try to reshape the data
try:
values_3d = data_isar_3d["Data"].to_numpy()
values_3d = values_3d.reshape((len(down_range), len(cross_range_az), len(cross_range_el)))
except Exception: # pragma: no cover
# fallback: fill with zeros
values_3d = np.zeros((len(down_range), len(cross_range_az), len(cross_range_el)))
if plane_cut.casefold() == "xy":
idx_fixed_range = (np.abs(cross_range_el - plane_offset)).argmin()
values = values_3d[:, :, idx_fixed_range]
y = down_range
x = cross_range_az
xlabel = "Down-Range (m)"
ylabel = "Cross Range (az) (m)"
elif plane_cut.casefold() == "xz":
idx_fixed_range = (np.abs(cross_range_az - plane_offset)).argmin()
values = values_3d[:, idx_fixed_range, :]
y = down_range
x = cross_range_el
xlabel = "Down-Range (m)"
ylabel = "Cross Range (el) (m)"
elif plane_cut.casefold() == "yz":
idx_fixed_range = (np.abs(down_range - plane_offset)).argmin()
values = values_3d[idx_fixed_range, :, :]
y = cross_range_az
x = cross_range_el
xlabel = "Cross Range (az) (m)"
ylabel = "Cross Range (el) (m)"
else:
raise ValueError("Invalid plane cut. Choose 'xy', 'xz', or 'yz'.")
x_, y_ = np.meshgrid(x, y) # important, do not use ij here
plot_data = [values, x_, y_]
new = ReportPlotter()
new.size = size
new.show_legend = False
new.title = title
props = {
"x_label": xlabel,
"y_label": ylabel,
}
new.add_trace(plot_data, 0, props)
_ = new.plot_pcolor(trace=0, snapshot_path=output_file, show=show, figure=figure)
return new
[docs]
@pyaedt_function_handler()
@graphics_required
def plot_scene(self, show=True):
"""Plot the 3D scene including models, annotations, and results.
This method visualizes the 3D scene by rendering the mesh objects under the "model",
"annotations", and "results" categories stored in `self.all_scene_actors`. The meshes
are rendered using a default PyVista plotter.
Parameters
----------
show : bool, default: ``True``
Whether to immediately display the plot using
the ``plotter.show()`` method. If ``False``, the
``plotter`` object is returned for further
customization before rendering.
Returns
-------
pyvista.Plotter or None
Returns the ``Plotter`` object if ``show`` is set to ``False``. If ``show`` is set to ``True``,
the plot is displayed and no value is returned.
"""
pv_backend = PyVistaBackend(allow_picking=True, plot_picked_names=True)
plotter = Plotter(backend=pv_backend)
if self.show_geometry:
for geo in self.all_scene_actors["model"].values():
self.__add_mesh(geo, plotter, "model")
for annotations in self.all_scene_actors["annotations"]:
for annotation in self.all_scene_actors["annotations"][annotations].values():
if annotation.custom_object.show:
self.__add_mesh(annotation, plotter, "annotations")
for all_scene_results in self.all_scene_actors["results"]:
for result_actor in self.all_scene_actors["results"][all_scene_results].values():
self.__add_mesh(result_actor, plotter, "results")
plotter.backend.scene.show_grid(
location="all",
xtitle=f"X Axis ({self.model_units})",
ytitle=f"Y Axis ({self.model_units})",
ztitle=f"Z Axis ({self.model_units})",
)
if show: # pragma: no cover
plotter.show()
else:
return plotter
[docs]
@pyaedt_function_handler()
@graphics_required
def add_rcs(
self,
color_bar="jet",
):
"""Add a 3D RCS representation to the current scene.
This function normalizes and visualizes RCS data on a spherical coordinate grid
(theta, phi), mapping it to 3D Cartesian coordinates (x, y, z). The RCS values are
color-mapped and added as a mesh to the current scene actors.
Parameters
----------
color_bar : str, default: ``"jet"``
Color mapping to apply to the RCS data. It can be a
color (such as ``"blue"`` or ``"green"``) or a
colormap (such as ``"jet"`` or ``"viridis"``.
"""
data = self.rcs_data.rcs_active_frequency
new_data = self.stretch_data(
data, scaling_factor=max(self.extents[5], 0.0) - min(self.extents[4], 0.0), offset=0.0
)
rcs = new_data["Data"]
rcs_min = np.min(rcs)
rcs_renorm = rcs + np.abs(rcs_min) if rcs_min else rcs
rcs_renorm = rcs_renorm.to_numpy()
# Negative values are not valid, this is cleaning numerical issues
rcs_renorm[rcs_renorm < 0] = 0.0
theta = np.deg2rad(data["IWaveTheta"])
phi = np.deg2rad(data["IWavePhi"])
unique_phi, unique_theta = np.unique(phi), np.unique(theta)
phi_grid, theta_grid = np.meshgrid(unique_phi, unique_theta)
r = np.full(phi_grid.shape, np.nan)
idx_theta = np.digitize(theta, unique_theta) - 1
idx_phi = np.digitize(phi, unique_phi) - 1
r[idx_theta, idx_phi] = rcs_renorm
x = r * np.sin(theta_grid) * np.cos(phi_grid)
y = r * np.sin(theta_grid) * np.sin(phi_grid)
z = r * np.cos(theta_grid)
actor = pv.StructuredGrid(x, y, z)
actor.point_data["values"] = data["Data"]
all_results_actors = list(self.all_scene_actors["results"].keys())
if "rcs" not in all_results_actors:
self.all_scene_actors["results"]["rcs"] = {}
index = 0
while f"rcs_{index}" in self.all_scene_actors["results"]["rcs"]:
index += 1
rcs_name = f"rcs_{index}"
rcs_object = SceneMeshObject()
rcs_object.name = rcs_name
rcs_object.line_width = 1.0
scalar_dict = dict(color="#000000", title=f"RCS {index} [{self.rcs_data.data_conversion_function}]")
rcs_object.scalar_dict = scalar_dict
if any(color_bar in x for x in ["blue", "green", "black", "red"]):
rcs_object.color = color_bar
else:
rcs_object.color_map = color_bar
rcs_object.mesh = actor
rcs_mesh = MeshObjectPlot(rcs_object, rcs_object.get_mesh())
self.all_scene_actors["results"]["rcs"][rcs_name] = rcs_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_range_profile_settings(
self,
size_range=10.0,
range_resolution=0.1,
tick_color="#000000",
line_color="#ff0000",
disc_color="ff0000",
cone_color="#00ff00",
):
"""Add a 3D range profile setting representation to the current scene.
This function visualizes a 3D range profile with a
main line representing the range axis and tick marks
indicating distance intervals. The profile includes
visual elements like a disc at the far end and a
cone at the starting point. These elements help to
display a reference range profile in the 3D scene.
Parameters
----------
size_range : float, default: ``10.0``
Total size of the range in meters. It determines the length of the range profile.
range_resolution : float, default: ``0.1``
Resolution of the range in meters, representing the distance between each tick mark along the range profile.
tick_color : str, default: ``"#000000"``
Color of the tick marks along the range profile. The ``"#000000"`` default is black.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` is red.
disc_color : str, default: ``"#ff0000"``
Color of the disc. The ``"#ff0000"`` is red.
cone_color : str, default: ``"#00ff00"``
Color of the cone. The d``"#00ff00"`` is green.
"""
size_range = unit_converter(
size_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
range_resolution = unit_converter(
range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
# Compute parameters
range_max = size_range - range_resolution
range_num = int(np.round(size_range / range_resolution))
distance_range = np.linspace(0, range_max, range_num)
distance_range -= distance_range[range_num // 2]
range_first = -distance_range[0]
range_last = -distance_range[-1]
num_ticks = int(size_range / range_resolution)
# Using 5% of total range length
tick_length = size_range * 0.05
if "range_profile" not in self.all_scene_actors["annotations"]:
self.all_scene_actors["annotations"]["range_profile"] = {}
# TODO: Do we want to support non-centered Range profile?
center = np.array([0.0, 0.0, 0.0])
# Main red line
name = "main_line"
main_line_az_mesh = self._create_line(
pointa=(range_first + center[0], self.extents[2] * 10 + center[1], center[2]),
pointb=(range_last + center[0], self.extents[2] * 10 + center[1], center[2]),
name=name,
color=line_color,
)
self.all_scene_actors["annotations"]["range_profile"][name] = main_line_az_mesh
# Ticks
tick_lines = pv.PolyData()
for tick in range(num_ticks + 1): # create line with tick marks
if tick % 1 == 0: # only do every nth tick
tick_pos_start = (
range_first - range_resolution * tick + center[0],
self.extents[2] * 10 + center[1],
center[2],
)
tick_pos_end = (
range_first - range_resolution * tick + center[0],
self.extents[2] * 10 + tick_length + center[1],
center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks"
tick_lines_object = SceneMeshObject()
tick_lines_object.name = annotation_name
tick_lines_object.color = tick_color
tick_lines_object.line_width = 2
tick_lines_object.mesh = tick_lines
tick_lines_mesh = MeshObjectPlot(tick_lines_object, tick_lines_object.get_mesh())
self.all_scene_actors["annotations"]["range_profile"][annotation_name] = tick_lines_mesh
start_geo = pv.Disc(
center=(range_last + center[0], self.extents[2] * 10 + center[1], center[2]),
outer=tick_length,
inner=0,
normal=(-1, 0, 0),
c_res=12,
)
annotation_name = "disc"
start_geo_object = SceneMeshObject()
start_geo_object.name = annotation_name
start_geo_object.color = disc_color
start_geo_object.line_width = 5
start_geo_object.mesh = start_geo
disc_geo_mesh = MeshObjectPlot(start_geo_object, start_geo_object.get_mesh())
self.all_scene_actors["annotations"]["range_profile"][annotation_name] = disc_geo_mesh
name = "cone"
cone_center = (range_first + center[0], self.extents[2] * 10 + center[1], center[2])
end_geo_mesh = self._create_cone(
center=cone_center,
direction=(-1, 0, 0),
radius=tick_length,
height=tick_length * 2,
resolution=12,
name=name,
color=cone_color,
)
self.all_scene_actors["annotations"]["range_profile"][name] = end_geo_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_waterfall_settings(
self, aspect_ang_phi=360.0, phi_num=10, tick_color="#000000", line_color="#ff0000", cone_color="#00ff00"
):
"""
Add a 3D waterfall setting representation to the current scene.
This function visualizes a 3D waterfall pattern that
represents angular data across a circular arc in a
spherical coordinate system. The arc covers an angular
extent defined by the ``aspect_ang_phi`` parameter, with
optional tick marks along the arc. A cone is added to
indicate the endpoint of the angular sweep.
Parameters
----------
aspect_ang_phi : float, default: ``360.0``
Angular extent of the arc in degrees. It defines
the total angle over which the circular arc spans.
The ``360.0`` default is a full circle.
phi_num : int, default: ``10``
Number of tick marks to place along the arc.
tick_color : str, default: ``"#000000"``
Color of the tick marks. The ``"#000000"`` default is black.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
cone_color : str, default: ``"#00ff00"``
Color of the cone. The ``"#00ff00"`` default is green.
"""
radius_max = self.radius
# TODO: Do we want to support non-centered waterfall?
center = np.array([0.0, 0.0, 0.0])
angle = aspect_ang_phi - 1
if "waterfall" not in self.all_scene_actors["annotations"]:
self.all_scene_actors["annotations"]["waterfall"] = {}
# Circular Arc
x_start = center[0] + radius_max
y_start = center[1]
start_point = (x_start, y_start, center[2])
end_point = [
center[0] + radius_max * np.cos(np.deg2rad(angle)),
center[1] + radius_max * np.sin(np.deg2rad(angle)),
center[2],
]
negative = False
if angle >= 180:
negative = True
name = "arc"
arc_mesh = self._create_arc(
pointa=start_point,
pointb=end_point,
center=center,
resolution=100,
negative=negative,
name="arc",
color=line_color,
)
self.all_scene_actors["annotations"]["waterfall"][name] = arc_mesh
# Ticks
tick_spacing_deg = int(aspect_ang_phi) / phi_num
if tick_spacing_deg >= 1: # gets too cluttered if too small of spacing
tick_lines = pv.PolyData()
for tick in range(phi_num + 1): # create line with tick marks
if tick % 1 == 0: # only do every nth tick
x_start = center[0] + radius_max * 0.95 * np.cos(np.deg2rad(tick * tick_spacing_deg))
y_start = center[1] + radius_max * 0.95 * np.sin(np.deg2rad(tick * tick_spacing_deg))
x_stop = center[0] + radius_max * 1.05 * np.cos(np.deg2rad(tick * tick_spacing_deg))
y_stop = center[1] + radius_max * 1.05 * np.sin(np.deg2rad(tick * tick_spacing_deg))
tick_pos_start = (x_start, y_start, center[2])
tick_pos_end = (x_stop, y_stop, center[2])
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks"
tick_lines_object = SceneMeshObject()
tick_lines_object.name = annotation_name
tick_lines_object.color = tick_color
tick_lines_object.line_width = 2
tick_lines_object.mesh = tick_lines
tick_lines_mesh = MeshObjectPlot(tick_lines_object, tick_lines_object.get_mesh())
self.all_scene_actors["annotations"]["waterfall"][annotation_name] = tick_lines_mesh
end_point = [
center[0] + radius_max * np.cos(np.deg2rad(angle)),
center[1] + radius_max * np.sin(np.deg2rad(angle)),
center[2],
]
end_point_plus_one = [
center[0] + radius_max * np.cos(np.deg2rad(aspect_ang_phi)),
center[1] + radius_max * np.sin(np.deg2rad(aspect_ang_phi)),
center[2],
]
direction = np.array(end_point_plus_one) - np.array(end_point)
direction_mag = np.linalg.norm(direction)
name = "cone"
end_geo_mesh = self._create_cone(
center=end_point,
direction=direction,
radius=direction_mag * 2,
height=direction_mag * 4,
resolution=12,
name=name,
color=cone_color,
)
self.all_scene_actors["annotations"]["waterfall"][name] = end_geo_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_isar_2d_settings(
self,
size_range=10.0,
range_resolution=0.1,
size_cross_range=10.0,
cross_range_resolution=0.1,
tick_color="#000000",
line_color="#ff0000",
):
"""
Add a preview frame of 2D ISAR visualization to the current 3D scene.
Parameters
----------
size_range : float, default: ``10.0``
Total size of the range axis in meters. This sets the overall length of the range axis.
range_resolution : float, default: ``0.1``
Resolution of the range axis in meters, specifying the
spacing between each tick mark.
size_cross_range : float, default: ``10.0``
Total size of the cross-range axis in meters. This sets the width of the cross-range axis.
cross_range_resolution : float, default: ``0.1``
Resolution of the cross-range axis in meters, specifying
the spacing between each tick mark along the azimuth axis.
tick_color : str, default: ``"#000000"``
Color of the tick marks along both the range and
cross-range axes. The "#000000" default is black.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
size_range = unit_converter(
size_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
range_resolution = unit_converter(
range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
size_cross_range = unit_converter(
size_cross_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
cross_range_resolution = unit_converter(
cross_range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
center = np.array([0.0, 0.0, 0.0])
range_ticks, range_frame, num_ticks = self._compute_axis_parameters(size_range, range_resolution, center[0])
range_ticks_az, range_frame_az, num_ticks_az = self._compute_axis_parameters(
size_cross_range, cross_range_resolution, center[1]
)
tick_length = size_range * 0.05
tick_length_az = size_cross_range * 0.05
if "isar_2d" not in self.all_scene_actors["annotations"]:
self.all_scene_actors["annotations"]["isar_2d"] = {}
# Create main lines
main_lines = self._create_main_lines(range_frame, range_frame_az, center, line_color)
self.all_scene_actors["annotations"]["isar_2d"].update(main_lines)
# Add range ticks
tick_lines = pv.PolyData()
for tick in range(0, num_ticks, num_ticks // 64 + 1): # create line with tick marks
tick_pos_start = (
range_ticks[tick] + center[0],
range_frame_az[0] + center[1],
center[2],
)
tick_pos_end = (
range_ticks[tick] + center[0],
range_frame_az[0] - tick_length + center[1],
center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks_range"
tick_lines_range_object = SceneMeshObject()
tick_lines_range_object.name = annotation_name
tick_lines_range_object.color = tick_color
tick_lines_range_object.line_width = 2
tick_lines_range_object.mesh = tick_lines
tick_lines_range_mesh = MeshObjectPlot(tick_lines_range_object, tick_lines_range_object.get_mesh())
self.all_scene_actors["annotations"]["isar_2d"][annotation_name] = tick_lines_range_mesh
# Add azimuth ticks
tick_lines = pv.PolyData()
for tick in range(0, num_ticks_az, num_ticks_az // 64 + 1): # create line with tick marks
tick_pos_start = (range_frame[-1] + center[0], range_ticks_az[tick] + center[1], center[2])
tick_pos_end = (
range_frame[-1] + tick_length_az + center[0],
range_ticks_az[tick] + center[1],
center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks_az"
tick_lines_az_object = SceneMeshObject()
tick_lines_az_object.name = annotation_name
tick_lines_az_object.color = tick_color
tick_lines_az_object.line_width = 2
tick_lines_az_object.mesh = tick_lines
tick_lines_az_mesh = MeshObjectPlot(tick_lines_az_object, tick_lines_az_object.get_mesh())
self.all_scene_actors["annotations"]["isar_2d"][annotation_name] = tick_lines_az_mesh
@staticmethod
def _compute_axis_parameters(size, resolution, center):
"""Compute axis parameters for ISAR settings."""
axis_max = size - resolution
axis_num = int(np.round(size / resolution))
axis_ticks = np.linspace(0, axis_max, axis_num)
axis_ticks -= (axis_ticks[-1] - axis_ticks[0]) / 2 + center
axis_frame = np.array([-size / 2, size / 2])
return axis_ticks, axis_frame, axis_num
def _create_main_lines(self, range_frame, range_frame_az, center, line_color):
"""Create main lines for ISAR 2D settings."""
main_lines = {}
main_lines["main_line"] = self._create_line(
pointa=(range_frame[0] + center[0], range_frame_az[0] + center[1], center[2]),
pointb=(range_frame[-1] + center[0], range_frame_az[0] + center[1], center[2]),
name="main_line",
color=line_color,
)
main_lines["main_line_opposite"] = self._create_line(
pointa=(range_frame[0] + center[0], range_frame_az[-1] + center[1], center[2]),
pointb=(range_frame[-1] + center[0], range_frame_az[-1] + center[1], center[2]),
name="main_line_opposite",
color=line_color,
)
main_lines["main_line_az"] = self._create_line(
pointa=(range_frame[0] + center[0], range_frame_az[0] + center[1], center[2]),
pointb=(range_frame[0] + center[0], range_frame_az[-1] + center[1], center[2]),
name="main_line_az",
color=line_color,
)
main_lines["main_line_az_opposite"] = self._create_line(
pointa=(range_frame[1] + center[0], range_frame_az[0] + center[1], center[2]),
pointb=(range_frame[1] + center[0], range_frame_az[1] + center[1], center[2]),
name="main_line_az_opposite",
color=line_color,
)
return main_lines
[docs]
@pyaedt_function_handler()
@graphics_required
def add_isar_3d_settings(
self,
size_range=10.0,
range_resolution=0.1,
size_cross_range=10.0,
cross_range_resolution=0.1,
size_elevation_range=10.0,
elevation_range_resolution=0.1,
tick_color="#000000",
line_color="#ff0000",
):
"""
Add a preview frame of 3D ISAR visualization to the current 3D scene.
Parameters
----------
size_range : float, default: ``10.0``
Total size of the range axis in meters. This sets the
overall length of the range axis.
range_resolution : float, default: ``0.1``
Resolution of the range axis in meters, specifying
the spacing between each tick mark.
size_cross_range : float, default: ``10.0``
Total size of the cross-range axis in meters. This
sets the width of the cross-range axis.
cross_range_resolution : float, default: ``0.1``
Resolution of the cross-range axis in meters, specifying
the spacing between each tick mark along the
azimuth axis.
size_elevation_range : float, default: ``10.0``
Total size of the elevation-range axis in meters. This
sets the width of the elevation-range axis.
elevation_range_resolution : float, default: ``0.1``
Resolution of the elevation-range axis in meters,
specifying the spacing between each tick mark along the elevation axis.
tick_color : str, default: ``"#000000"``
Color of the tick marks along both the range and
cross-range axes. The ``"#000000"`` default is black.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
size_range = unit_converter(
size_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
range_resolution = unit_converter(
range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
size_cross_range = unit_converter(
size_cross_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
cross_range_resolution = unit_converter(
cross_range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
size_elevation_range = unit_converter(
size_elevation_range, unit_system="Length", input_units="meter", output_units=self.model_units
)
elevation_range_resolution = unit_converter(
elevation_range_resolution, unit_system="Length", input_units="meter", output_units=self.model_units
)
center = np.array([0.0, 0.0, 0.0])
range_ticks, range_frame, num_ticks = self._compute_axis_parameters(size_range, range_resolution, center[0])
range_ticks_az, range_frame_az, num_ticks_az = self._compute_axis_parameters(
size_cross_range, cross_range_resolution, center[1]
)
range_ticks_el, range_frame_el, num_ticks_el = self._compute_axis_parameters(
size_elevation_range, elevation_range_resolution, center[2]
)
tick_length = size_range * 0.05
tick_length_az = size_cross_range * 0.05
tick_length_el = size_elevation_range * 0.05
if "isar_3d" not in self.all_scene_actors["annotations"]:
self.all_scene_actors["annotations"]["isar_3d"] = {}
# Create main lines
main_lines = self._create_main_lines(range_frame, range_frame_az, center, line_color)
self.all_scene_actors["annotations"]["isar_3d"].update(main_lines)
# Add elevation lines
name = "main_line_el"
main_line_el_mesh = self._create_line(
pointa=(range_frame[0] + center[0], range_frame_az[0] + center[1], range_frame_el[0] + center[2]),
pointb=(range_frame[0] + center[0], range_frame_az[0] + center[1], range_frame_el[-1] + center[2]),
name=name,
color=line_color,
)
self.all_scene_actors["annotations"]["isar_3d"][name] = main_line_el_mesh
# Add bounding box
box_bounds = (
range_frame[0] + center[0],
range_frame[-1] + center[0],
range_frame_az[0] + center[1],
range_frame_az[-1] + center[1],
range_frame_el[0] + center[2],
range_frame_el[-1] + center[2],
)
name = "box"
box = pv.Box(box_bounds)
box_object = SceneMeshObject()
box_object.name = name
box_object.color = line_color
box_object.opacity = 0.05
box_object.line_width = 5
box_object.show_edges = True
box_object.edge_color = line_color
box_object.mesh = box
box_mesh = MeshObjectPlot(box_object, box_object.get_mesh())
self.all_scene_actors["annotations"]["isar_3d"][name] = box_mesh
# Add range ticks
if num_ticks <= 256:
tick_lines = pv.PolyData()
for tick in range(1, num_ticks, 2):
tick_pos_start = (
range_ticks[tick] + center[0],
range_frame_az[0] + center[1],
center[2],
)
tick_pos_end = (
range_ticks[tick] + center[0],
range_frame_az[0] - tick_length + center[1],
center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks_range"
tick_lines_range_object = SceneMeshObject()
tick_lines_range_object.name = annotation_name
tick_lines_range_object.color = tick_color
tick_lines_range_object.line_width = 2
tick_lines_range_object.mesh = tick_lines
tick_lines_range_mesh = MeshObjectPlot(tick_lines_range_object, tick_lines_range_object.get_mesh())
self.all_scene_actors["annotations"]["isar_3d"][annotation_name] = tick_lines_range_mesh
# Add azimuth ticks
if num_ticks_az <= 256:
tick_lines = pv.PolyData()
for tick in range(1, num_ticks_az - 1, 2):
tick_pos_start = (range_frame[0] + center[0], range_ticks_az[tick] + center[1], center[2])
tick_pos_end = (
range_frame[0] + tick_length_az + center[0],
range_ticks_az[tick] + center[1],
center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks_az"
tick_lines_az_object = SceneMeshObject()
tick_lines_az_object.name = annotation_name
tick_lines_az_object.color = tick_color
tick_lines_az_object.line_width = 2
tick_lines_az_object.mesh = tick_lines
tick_lines_az_mesh = MeshObjectPlot(tick_lines_az_object, tick_lines_az_object.get_mesh())
self.all_scene_actors["annotations"]["isar_3d"][annotation_name] = tick_lines_az_mesh
# Add elevation ticks
if num_ticks_el <= 256:
tick_lines = pv.PolyData()
for tick in range(1, num_ticks_el - 1, 2):
tick_pos_start = (
range_frame[0] + center[0],
range_frame_az[0] + center[1],
range_ticks_el[tick] + center[2],
)
tick_pos_end = (
range_frame[0] + tick_length_el + center[0],
range_frame_az[0] + center[1],
range_ticks_el[tick] + center[2],
)
tick_lines += pv.Line(pointa=tick_pos_start, pointb=tick_pos_end)
annotation_name = "ticks_el"
tick_lines_el_object = SceneMeshObject()
tick_lines_el_object.name = annotation_name
tick_lines_el_object.color = tick_color
tick_lines_el_object.line_width = 2
tick_lines_el_object.mesh = tick_lines
tick_lines_el_mesh = MeshObjectPlot(tick_lines_el_object, tick_lines_el_object.get_mesh())
self.all_scene_actors["annotations"]["isar_3d"][annotation_name] = tick_lines_el_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_range_profile(
self,
plot_type="Line",
color_bar="jet",
):
"""Add the 3D range profile.
This function visualizes a 3D range profile, which represents the RCS data
as a function of range. The profile is rendered as a plot in the 3D scene using spherical coordinates
(phi, theta) mapped into Cartesian coordinates (x, y). The range profile can be visualized in different
plot types such as line plots or color-mapped surfaces.
Parameters
----------
plot_type : str, default: ``"Line"``
The type of plot to create for the range profile.
Options are ``"Extruded"``, ``"Line"``, `"Plane H"``,
``"Plane V"``, `"Projection"``, ``"Ribbon"``, and ``"Rotated"``.
color_bar : str, default: ``"jet"``
Color mapping to apply to the RCS data. It can be
a color (such as ``"blue"`` or ``"green"``) or a
colormap (such as ``"jet"`` or ``"viridis"``).
"""
data_range_profile = self.rcs_data.range_profile
data_scaled = self.stretch_data(
data_range_profile["Data"],
scaling_factor=max(self.extents[5], 0.0) - min(self.extents[4], 0.0),
offset=0.0,
)
ranges = data_range_profile["Range"].to_numpy()
rotation = {}
rotation["azimuth"] = -self.rcs_data.incident_wave_phi
rotation["elevation"] = 90 - self.rcs_data.incident_wave_theta
rotation["twist"] = 0
xpos = -ranges
ypos = np.zeros_like(xpos)
if plot_type.casefold() in ["line", "ribbon", "rotated", "extruded"]:
zpos = data_scaled.to_numpy()
else:
zpos = np.zeros_like(xpos)
x, y, z = self.rotate_point(xpos, ypos, zpos, rotation)
cpos = data_range_profile["Data"].to_numpy()
vmin, vmax = np.nanmin(cpos), np.nanmax(cpos)
xvec, yvec, zvec = self.rotate_point(np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1]), rotation)
axis_vec = np.array([[xvec[0], yvec[0], zvec[0]], [xvec[1], yvec[1], zvec[1]], [xvec[2], yvec[2], zvec[2]]])
actor = self.__get_pyvista_range_profile_actor(
x,
y,
z,
cpos,
axis_vec,
plot_type=plot_type,
scene_actors=self.all_scene_actors["model"],
data_conversion_function=self.rcs_data.data_conversion_function,
extents=self.extents,
)
all_results_actors = list(self.all_scene_actors["results"].keys())
if "range_profile" not in all_results_actors:
self.all_scene_actors["results"]["range_profile"] = {}
index = 0
while f"range_profile_{index}" in self.all_scene_actors["results"]["range_profile"]:
index += 1
range_profile_name = f"range_profile_{index}"
range_profile_object = SceneMeshObject()
range_profile_object.name = range_profile_name
range_profile_object.line_width = 1.0
scalar_dict = dict(color="#000000", title=f"Range Profile {index} [{self.rcs_data.data_conversion_function}]")
range_profile_object.scalar_dict = scalar_dict
if any(color_bar in x for x in ["blue", "green", "black", "red"]):
range_profile_object.color = color_bar
else:
range_profile_object.color_map = color_bar
range_profile_object.mesh = actor
range_profile_object.clim = (vmin, vmax)
range_profile_object.default_clim = (vmin, vmax)
range_profile_mesh = MeshObjectPlot(range_profile_object, range_profile_object.get_mesh())
self.all_scene_actors["results"]["range_profile"][range_profile_name] = range_profile_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_waterfall(
self,
color_bar="jet",
):
"""Add the 3D waterfall."""
data_waterfall = self.rcs_data.waterfall
rotation = {}
if self.rcs_data.aspect_range == "Horizontal":
rotation["azimuth"] = 0
rotation["elevation"] = 90 - self.rcs_data.incident_wave_theta
rotation["twist"] = 0
angles = np.unique(data_waterfall["IWavePhi"]) # negative azimuth
elif self.rcs_data.aspect_range == "Vertical":
rotation["azimuth"] = -self.rcs_data.incident_wave_phi
rotation["elevation"] = 0
rotation["twist"] = 90
angles = np.unique(data_waterfall["IWaveTheta"]) - 90 # negative elevation
else:
raise ValueError("Plane selection is not supported. Please select 'Horizontal' or 'Vertical'.")
ranges = np.unique(data_waterfall["Range"])
values = data_waterfall["Data"].to_numpy()
angles = np.deg2rad(angles)
ra, ph = np.meshgrid(ranges, angles)
x = -ra.T * np.cos(ph.T)
y = -ra.T * np.sin(ph.T)
z = np.zeros_like(x)
# Rotate the points based on the specified rotation
x, y, z = self.rotate_point(x, y, z, rotation)
actor = pv.StructuredGrid(x, y, z)
actor.point_data["values"] = values
all_results_actors = list(self.all_scene_actors["results"].keys())
if "waterfall" not in all_results_actors:
self.all_scene_actors["results"]["waterfall"] = {}
index = 0
while f"waterfall_{index}" in self.all_scene_actors["results"]["waterfall"]:
index += 1
waterfall_name = f"waterfall_{index}"
waterfall_object = SceneMeshObject()
waterfall_object.name = waterfall_name
scalar_dict = dict(color="#000000", title=f"Waterfall {index} [{self.rcs_data.data_conversion_function}]")
waterfall_object.scalar_dict = scalar_dict
waterfall_object.color_map = color_bar
waterfall_object.mesh = actor
rcs_mesh = MeshObjectPlot(waterfall_object, waterfall_object.get_mesh())
self.all_scene_actors["results"]["waterfall"][waterfall_name] = rcs_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_isar_2d(
self,
plot_type="plane",
color_bar="jet",
):
"""Add the ISAR 2D.
Parameters
----------
plot_type : str, default: ``"plane"``
Type of plot to create for the range profile. Options
are ``"plane"``, `"projection"`` and ``"relief"``.
color_bar : str, default: ``"jet"``
Color mapping to apply to the RCS data. It can be a
color (such as ``"blue"`` or ``"green"``) or a colormap
(such as ``"jet"`` or ``"viridis"``).
"""
rotation = {}
if self.rcs_data.aspect_range == "Horizontal":
rotation["azimuth"] = 0
rotation["elevation"] = 90 - self.rcs_data.incident_wave_theta
rotation["twist"] = 0
elif self.rcs_data.aspect_range == "Vertical":
rotation["azimuth"] = -self.rcs_data.incident_wave_phi
rotation["elevation"] = 0
rotation["twist"] = 90
else:
raise ValueError("Plane selection is not supported. Please select 'Horizontal' or 'Vertical'.")
data_isar_2d = self.rcs_data.isar_2d
down_range = data_isar_2d["Down-range"].unique()
cross_range = data_isar_2d["Cross-range"].unique()
values_2d = data_isar_2d["Data"].to_numpy().reshape((len(down_range), len(cross_range)))
if plot_type.casefold() in ["relief", "plane"]:
# meshgrid must have one more pixel. In the other words, number of meshgrid points = number of values + 1
# mesh idx 1 2 3 4 5
# | * | * | * | * |
# value idx 1 2 3 4
dx = down_range[1] - down_range[0]
down_range_grid = np.linspace(down_range[0] - dx / 2, down_range[-1] + dx / 2, num=len(down_range) + 1)
dy = cross_range[1] - cross_range[0]
cross_range_grid = np.linspace(cross_range[0] - dy / 2, cross_range[-1] + dy / 2, num=len(cross_range) + 1)
x, y = np.meshgrid(down_range_grid[::-1], cross_range_grid, indexing="ij")
z = np.zeros_like(x)
if plot_type.casefold() == "relief":
m = 2.0
b = -1.0
# z = (values_2d - values_2d.min()) / (values_2d.max() - values_2d.min()) * m + b
f_z = RegularGridInterpolator(
(-down_range, cross_range), values_2d, fill_value=None, method="linear", bounds_error=False
)
z = f_z((x, y))
z = (z - z.min()) / (z.max() - z.min()) * m + b
else:
# Projection
x, y = np.meshgrid(down_range[::-1], cross_range, indexing="ij")
z = np.zeros_like(x)
# Rotate the points based on the specified rotation
x, y, z = self.rotate_point(x, y, z, rotation)
vmin, vmax = np.nanmin(values_2d), np.nanmax(values_2d)
if plot_type.casefold() in ["relief", "plane"]:
actor = pv.StructuredGrid(x, y, z)
actor["values"] = values_2d.ravel(order="F")
else: # Projection
scene_actors = self.all_scene_actors["model"]
if scene_actors is None: # pragma: no cover
return None
actor = pv.PolyData()
# Prepare grid geometry for fast nearest neighbor search
# x, y, z are meshgrid arrays (down_range, cross_range, 0), shape (N, M)
# The grid is in the (x, y) plane, z=0
# Grid origin is the (x[0,0], y[0,0], z[0,0])
cpos = values_2d.flatten()
isar_points = np.column_stack((x.flatten(), y.flatten(), z.flatten()))
grid_origin = isar_points[0, :]
# Basis vectors: u = (x[1,0] - x[0,0], y[1,0] - y[0,0], z[1,0] - z[0,0])
# v = (x[0,1] - x[0,0], y[0,1] - y[0,0], z[0,1] - z[0,0])
if x.shape[0] > 1:
u = np.array([x[1, 0] - x[0, 0], y[1, 0] - y[0, 0], z[1, 0] - z[0, 0]])
else: # pragma: no cover
u = np.array([1.0, 0.0, 0.0])
if x.shape[1] > 1:
v = np.array([x[0, 1] - x[0, 0], y[0, 1] - y[0, 0], z[0, 1] - z[0, 0]])
else: # pragma: no cover
v = np.array([0.0, 1.0, 0.0])
grid_shape = x.shape
for model_actor in scene_actors.values():
mesh = model_actor.custom_object.get_mesh()
target_points = mesh.points
all_indices = self.__find_nearest_neighbors(
isar_points, target_points, grid_origin=grid_origin, u=u, v=v, grid_shape=grid_shape
)
mag_for_color = np.ndarray.flatten(cpos[all_indices])
if mesh.__class__.__name__ != "PolyData": # pragma: no cover
mesh_triangulated = mesh.triangulate()
model_actor.custom_object.mesh = pv.PolyData(mesh_triangulated.points, mesh_triangulated.cells)
else:
model_actor.custom_object.mesh.clear_data()
model_actor.custom_object.mesh[self.rcs_data.data_conversion_function] = mag_for_color
actor += model_actor.custom_object.mesh
all_results_actors = list(self.all_scene_actors["results"].keys())
if "isar_2d" not in all_results_actors:
self.all_scene_actors["results"]["isar_2d"] = {}
index = 0
while f"isar_2d_{index}" in self.all_scene_actors["results"]["isar_2d"]:
index += 1
isar_name = f"isar_2d_{index}"
isar_object = SceneMeshObject()
isar_object.name = isar_name
scalar_dict = dict(color="#000000", title=f"ISAR 2D {index} [{self.rcs_data.data_conversion_function}]")
isar_object.scalar_dict = scalar_dict
isar_object.color_map = color_bar
isar_object.mesh = actor
isar_object.clim = (vmin, vmax)
isar_object.default_clim = (vmin, vmax)
isar_object.plot_type = plot_type.lower()
rcs_mesh = MeshObjectPlot(isar_object, isar_object.get_mesh())
self.all_scene_actors["results"]["isar_2d"][isar_name] = rcs_mesh
[docs]
@pyaedt_function_handler()
@graphics_required
def add_isar_3d(
self, plot_type: str = "iso-surface", color_bar: str = "jet", plane_cut: str = None, plane_offset: float = None
) -> None:
"""Add the ISAR 3D plot to the 3D scene.
Parameters
----------
plot_type : str, default: ``"iso-surface"``
Type of plot to create for the range profile. Options are `"iso-surface"`, `"plane cut"`, `"point cloud"`,
and `"projection"`.
color_bar : str, default: ``"jet"``
Color mapping to apply to the RCS data. It can be a
color (such as ``"blue"`` or ``"green"``) or a colormap
(such as ``"jet"``or ``"viridis"``).
plane_cut : str, default: ``None``
Type of plane cut to create for 3D ISAR. Options are
``"xy"``, ``"xz"``, and ``"yz"``. The ``None`` default
is for other plot types.
plane_offset : float, default: ``None``
Offset value for the plane cut, which is used to
specify the position of the plane in the 3D ISAR
data. The ``None`` default is for other plot types.
"""
data_isar_3d = self.rcs_data.isar_3d
down_range = data_isar_3d["Down-range"].unique()
cross_range_az = data_isar_3d["Cross-range-az"].unique()
cross_range_el = data_isar_3d["Cross-range-el"].unique()
dx = down_range[1] - down_range[0] if len(down_range) > 1 else 0
dy = cross_range_az[1] - cross_range_az[0] if len(cross_range_az) > 1 else 0
dz = cross_range_el[1] - cross_range_el[0] if len(cross_range_el) > 1 else 0
down_range_grid = np.linspace(down_range[0] - dx / 2, down_range[-1] + dx / 2, num=len(down_range) + 1)
cross_range_az_grid = np.linspace(
cross_range_az[0] - dy / 2, cross_range_az[-1] + dy / 2, num=len(cross_range_az) + 1
)
cross_range_el_grid = np.linspace(
cross_range_el[0] - dz / 2, cross_range_el[-1] + dz / 2, num=len(cross_range_el) + 1
)
# The above pivot may need to be adjusted depending on the data layout
# For now, let's try to reshape the data
# instead of modifying the direction of the down_range grid, we just flip the order of the data
# isar 3d. down_range now means "x" axis.
try:
values_3d = data_isar_3d["Data"].to_numpy()
values_3d = values_3d.reshape((len(down_range), len(cross_range_az), len(cross_range_el)))
values_3d = values_3d[::-1, :, :] # reverse down_range to match the original data orientation
except Exception: # pragma: no cover
# fallback: fill with zeros
values_3d = np.zeros((len(down_range), len(cross_range_az), len(cross_range_el)))
all_results_actors = list(self.all_scene_actors["results"].keys())
if "isar_3d" not in all_results_actors:
self.all_scene_actors["results"]["isar_3d"] = {}
index = 0
while f"isar_3d_{index}" in self.all_scene_actors["results"]["isar_3d"]:
index += 1
isar_name = f"isar_3d_{index}"
vmin, vmax = np.nanmin(values_3d), np.nanmax(values_3d)
if plot_type.casefold() != "plane cut":
# Implementation for point cloud, iso-surface, and projection plots
if plot_type.casefold() == "point cloud":
shape = (len(down_range_grid), len(cross_range_az_grid), len(cross_range_el_grid))
spacing = [dx, dy, dz]
origin = (np.min(down_range_grid), np.min(cross_range_az_grid), np.min(cross_range_el_grid))
else:
# iso-surface and projection require the center of the points, not the grid
shape = (len(down_range), len(cross_range_az), len(cross_range_el))
spacing = [dx, dy, dz]
origin = (np.min(down_range), np.min(cross_range_az), np.min(cross_range_el))
else:
# implementation for plane cut plot
if plane_cut is None:
raise ValueError("Please provide a plane cut type.")
if plane_offset is None:
raise ValueError("Please provide a plane offset value.")
if plane_cut.casefold() == "xy":
idx_fixed_range = (np.abs(cross_range_el - plane_offset)).argmin()
values_3d = values_3d[:, :, idx_fixed_range]
shape = (len(down_range_grid), len(cross_range_az_grid), 1)
spacing = [dx, dy, 0]
origin = (np.min(down_range_grid), np.min(cross_range_az_grid), cross_range_el[idx_fixed_range])
elif plane_cut.casefold() == "xz":
idx_fixed_range = (np.abs(cross_range_az - plane_offset)).argmin()
values_3d = values_3d[:, idx_fixed_range, :]
shape = (len(down_range_grid), 1, len(cross_range_el_grid))
spacing = [dx, 0, dz]
origin = (np.min(down_range_grid), cross_range_az[idx_fixed_range], np.min(cross_range_el_grid))
elif plane_cut.casefold() == "yz":
idx_fixed_range = (np.abs(down_range + plane_offset)).argmin() # down_range is negative
values_3d = values_3d[idx_fixed_range, :, :]
shape = (1, len(cross_range_az_grid), len(cross_range_el_grid))
spacing = [0, dy, dz]
origin = (down_range[idx_fixed_range], np.min(cross_range_az_grid), np.min(cross_range_el_grid))
else:
raise ValueError("Invalid plane cut. Choose 'xy', 'xz', or 'yz'.")
if plot_type.casefold() != "projection":
actor = pv.ImageData(dimensions=shape, spacing=spacing, origin=origin)
actor["scalars"] = values_3d.ravel(order="F")
else:
# configuration of this actor will happen in the relevant code block
actor = pv.PolyData()
isar_object = SceneMeshObject()
isar_object.name = isar_name
scalar_dict = dict(color="#FFFFFF", title=f"ISAR 3D {index} [{self.rcs_data.data_conversion_function}]")
isar_object.scalar_dict = scalar_dict
isar_object.color_map = color_bar
if plot_type.casefold() == "point cloud":
isar_object.mesh = actor
isar_object.object_type = SceneMeshObjectType.VOLUME
isar_object.opacity = "linear" if "dB" not in self.rcs_data.data_conversion_function else "geom"
elif plot_type.casefold() == "plane cut":
isar_object.mesh = actor
isar_object.object_type = SceneMeshObjectType.MESH
isar_object.opacity = 1.0
elif plot_type.casefold() == "iso-surface":
contour_values = np.linspace(vmin, vmax, self.num_contours)
contours = actor.contour(isosurfaces=contour_values)
isar_object.mesh = contours
isar_object.object_type = SceneMeshObjectType.MESH
isar_object.opacity = "linear" if "dB" not in self.rcs_data.data_conversion_function else "geom"
else: # Projection
scene_actors = self.all_scene_actors["model"]
if scene_actors is None: # pragma: no cover
return None
# Prepare grid geometry for fast nearest neighbor search
# x, y, z are meshgrid arrays (down_range, cross_range, 0), shape (N, M)
# The grid is in the (x, y) plane, z=0
# Grid origin is the (x[0,0], y[0,0], z[0,0])
x, y, z = np.meshgrid(down_range, cross_range_az, cross_range_el, indexing="ij")
cpos = values_3d.flatten()
isar_points = np.column_stack((x.flatten(), y.flatten(), z.flatten()))
grid_origin = isar_points[0, :]
# Basis vectors: u = (x[1,0,0] - x[0,0,0], y[1,0,0] - y[0,0,0], z[1,0,0] - z[0,0,0])
# v = (x[0,1,0] - x[0,0,0], y[0,1,0] - y[0,0,0], z[0,1,0] - z[0,0,0])
# w = (x[0,0,1] - x[0,0,0], y[0,0,1] - y[0,0,0], z[0,0,1] - z[0,0,0])
if x.shape[0] > 1:
u = np.array([x[1, 0, 0] - x[0, 0, 0], y[1, 0, 0] - y[0, 0, 0], z[1, 0, 0] - z[0, 0, 0]])
else: # pragma: no cover
u = np.array([1.0, 0.0, 0.0])
if x.shape[1] > 1:
v = np.array([x[0, 1, 0] - x[0, 0, 0], y[0, 1, 0] - y[0, 0, 0], z[0, 1, 0] - z[0, 0, 0]])
else: # pragma: no cover
v = np.array([0.0, 1.0, 0.0])
if x.shape[2] > 1:
w = np.array([x[0, 0, 1] - x[0, 0, 0], y[0, 0, 1] - y[0, 0, 0], z[0, 0, 1] - z[0, 0, 0]])
else: # pragma: no cover
w = np.array([0.0, 0.0, 1.0])
grid_shape = x.shape
for model_actor in scene_actors.values():
mesh = model_actor.custom_object.get_mesh()
target_points = mesh.points
all_indices = self.__find_nearest_neighbors(
isar_points, target_points, grid_origin=grid_origin, u=u, v=v, w=w, grid_shape=grid_shape
)
mag_for_color = np.ndarray.flatten(cpos[all_indices])
if mesh.__class__.__name__ != "PolyData": # pragma: no cover
mesh_triangulated = mesh.triangulate()
model_actor.custom_object.mesh = pv.PolyData(mesh_triangulated.points, mesh_triangulated.cells)
else:
model_actor.custom_object.mesh.clear_data()
model_actor.custom_object.mesh[self.rcs_data.data_conversion_function] = mag_for_color
actor += model_actor.custom_object.mesh
isar_object.mesh = actor
isar_object.object_type = SceneMeshObjectType.MESH
isar_object.opacity = 1.0
isar_object.clim = (vmin, vmax)
isar_object.default_clim = (vmin, vmax)
isar_object.plot_type = plot_type.lower()
isar_object.default_mesh = actor
rcs_mesh = MeshObjectPlot(isar_object, isar_object.get_mesh())
self.all_scene_actors["results"]["isar_3d"][isar_name] = rcs_mesh
[docs]
def get_plane_cut_from_3d_isar(self, shape, trace_properties): # pragma: no cover
"""
Extract a 2D plane cut from 3D ISAR data based on the specified plane.
This method computes a slice of the 3D ISAR dataset by selecting the closest available
plane offset in the data. It returns the corresponding sliced data and the new shape
of the 2D array.
Parameters
----------
shape : tuple of int
The shape of the original 3D ISAR data array.
trace_properties : dict
Dictionary containing the following keys:
``"plane_cut"``: str.
Indicates the slicing plane. Must be one of ``"XY-Plane"``, ``"XZ-Plane"``, or ``"YZ-Plane"``.
``"plane_offset"``: float.
The target coordinate value at which to slice the selected plane.
Returns
-------
new_data : pandas.DataFrame
A 2D slice of the ISAR data at the nearest available offset.
new_shape : tuple of int
Shape of the resulting 2D array after slicing. One of the dimensions will be ``1``.
"""
plane_offset = trace_properties["plane_offset"]
if trace_properties["plane_cut"] == "XY-Plane":
slice_range_name = "Cross-range2"
new_shape = (shape[0], shape[1], 1)
elif trace_properties["plane_cut"] == "XZ-Plane":
slice_range_name = "Cross-range1"
new_shape = (shape[0], 1, shape[2])
else:
slice_range_name = "Down-range"
new_shape = (1, shape[1], shape[2])
slice_range = self.rcs_data.raw_data.index.get_level_values(slice_range_name)
idx_fixed_range = (np.abs(slice_range - plane_offset)).argmin()
trace_properties["plane_offset"] = slice_range[idx_fixed_range] # replace with actual value
new_data = self.rcs_data.raw_data.loc[slice_range == trace_properties["plane_offset"]]
return new_data, new_shape
[docs]
@pyaedt_function_handler()
def add_incident_rcs_settings(
self,
theta_span: float,
num_theta: int,
phi_span: float,
num_phi: int,
arrow_color: str = "#ff0000",
line_color: str = "#ff0000",
):
"""Add incident wave arrow setting for RCS scene.
This function visualizes the incident wave arrows for RCS settings.
Parameters
----------
theta_span : float
Incident theta angle in degrees.
num_theta : int
Number of theta points.
phi_span : float
Incident phi angle in degrees.
num_phi : int
Number of phi points.
arrow_color : str, default: ``"#ff0000"``
Color of the arrow. The ``"#ff0000"`` default is red.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
self._add_incident_settings(
scene_type="rcs",
theta_span=theta_span,
num_theta=num_theta,
phi_span=phi_span,
num_phi=num_phi,
arrow_color=arrow_color,
line_color=line_color,
)
[docs]
@pyaedt_function_handler()
def add_incident_range_profile_settings(self, arrow_color: str = "#ff0000"):
"""Add incident wave arrow setting for range profile scene.
This function visualizes the incident wave arrows for RCS settings.
Parameters
----------
arrow_color : str, default: ``"#ff0000"``
Color of the arrow. The ``"#ff0000"`` default is red.
"""
self._add_incident_settings(scene_type="range_profile", arrow_color=arrow_color)
[docs]
@pyaedt_function_handler()
def add_incident_waterfall_settings(
self, phi_span: float, num_phi: int, arrow_color: str = "#ff0000", line_color: str = "#ff0000"
):
"""Add incident wave arrow setting for waterfall scene.
This function visualizes the incident wave arrows for waterfall settings.
Parameters
----------
phi_span : float
Incident phi angle in degrees.
num_phi : int
Number of phi points.
arrow_color : str, default: ``"#ff0000"``
Color of the arrow. The ``"#ff0000"`` default is red.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
self._add_incident_settings(
scene_type="waterfall", phi_span=phi_span, num_phi=num_phi, arrow_color=arrow_color, line_color=line_color
)
[docs]
@pyaedt_function_handler()
def add_incident_isar_2d_settings(
self, phi_span: float, num_phi: int, arrow_color: str = "#ff0000", line_color: str = "#ff0000"
):
"""Add incident wave arrow setting for ISAR 2D scene.
This function visualizes the incident wave arrows for ISAR 2D settings.
Parameters
----------
phi_span : float
Incident phi angle in degrees.
num_phi : int
Number of phi points.
arrow_color : str, default: ``"#ff0000"``
Color of the arrow. The ``"#ff0000"`` default is red.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
self._add_incident_settings(
scene_type="isar_2d", phi_span=phi_span, num_phi=num_phi, arrow_color=arrow_color, line_color=line_color
)
[docs]
@pyaedt_function_handler()
def add_incident_isar_3d_settings(
self,
theta_span: float,
num_theta: int,
phi_span: float,
num_phi: int,
arrow_color: str = "#ff0000",
line_color: str = "#ff0000",
):
"""Add incident wave arrow setting for ISAR 3D scene.
This function visualizes the incident wave arrows for ISAR 3D settings.
Parameters
----------
theta_span : float
Incident theta angle in degrees.
num_theta : int
Number of theta points.
phi_span : float
Incident phi angle in degrees.
num_phi : int
Number of phi points.
arrow_color : str, default: ``"#ff0000"``
Color of the arrow. The ``"#ff0000"`` default is red.
line_color : str, default: ``"#ff0000"``
Color of the line. The ``"#ff0000"`` default is red.
"""
self._add_incident_settings(
scene_type="isar_3d",
theta_span=theta_span,
num_theta=num_theta,
phi_span=phi_span,
num_phi=num_phi,
arrow_color=arrow_color,
line_color=line_color,
)
[docs]
@pyaedt_function_handler()
def clear_scene(self, first_level=None, second_level=None, name=None):
"""
Clear actors from the visualization scene.
This method removes elements from the scene based on the provided hierarchy levels.
By default, it clears all annotations and results. The geometry (model) cannot be cleared using this method.
Parameters
----------
first_level : str, default: ``None``
Top-level key indicating what to clear. Options are
``"annotations"``, ``"model"``, and ``"results"``.
If no option is provided, both annotations and
results are cleared.
second_level : str, default: ``None``
Subcategory within the specified ``first_level``
key. If no value is provided, all data under the ``first_level`` key is cleared.
name : str, default: ``None``
Specific item to remove within the ``second_level`` subcategory. If no value is provided,
all items under the ``second_level`` subcategory
are removed.
Returns
-------
bool
Returns ``True`` successful, ``False`` otherwise.
"""
if not first_level:
self.all_scene_actors["annotations"] = {}
self.all_scene_actors["results"] = {}
elif first_level == "model":
self.__logger.warning("Model can not be cleared. Set 'show_geometry' to False.")
return False
elif first_level in ["annotations", "results"]:
if not second_level:
self.all_scene_actors[first_level] = {}
elif second_level in self.all_scene_actors[first_level].keys(): # pragma: no cover
if not name:
self.all_scene_actor[first_level][second_level] = {}
elif name in self.all_scene_actors[first_level][second_level].keys():
del self.all_scene_actors[first_level][second_level][name]
return True
@pyaedt_function_handler()
def _add_incident_settings(
self,
scene_type="RCS",
theta_span=0.0,
num_theta=101,
phi_span=0.0,
num_phi=101,
arrow_color="#ff0000",
line_color="#ff0000",
):
# Compute parameters
radius_max = max([abs(self.extents[0]), abs(self.extents[1]), abs(self.extents[2]), abs(self.extents[3])])
arrow_length = radius_max * 0.25
if "incident_wave" not in self.all_scene_actors["annotations"]:
self.all_scene_actors["annotations"]["incident_wave"] = {}
# Plot arrows
# if we end up wanting to plot points, we need to keep the whole span.
# Arrows are only plotted at the edges of the domain
theta = np.linspace(0, theta_span, num_theta)
theta -= theta[len(theta) // 2]
theta += 90
theta = np.deg2rad([theta[0], theta[-1]])
phi = np.linspace(0, phi_span, num_phi)
phi -= phi[len(phi) // 2]
phi = np.deg2rad([phi[0], phi[-1]])
corners = ((theta[0], phi[0]), (theta[0], phi[1]), (theta[1], phi[1]), (theta[1], phi[0]))
for i, (t, p) in enumerate(corners):
arrow_direction = -np.array([np.cos(p) * np.sin(t), np.sin(p) * np.sin(t), np.cos(t)])
arrow_start = -arrow_direction * (radius_max + arrow_length)
name = "arrow" + str(i)
arrow_mesh = self._create_arrow(
start=arrow_start, direction=arrow_direction, scale=arrow_length, name=name, color=arrow_color
)
self.all_scene_actors["annotations"]["incident_wave"][name] = arrow_mesh
corner1 = arrow_start
n_t, n_p = corners[(i + 1) % 4]
n_arrow_direction = -np.array([np.cos(n_p) * np.sin(n_t), np.sin(n_p) * np.sin(n_t), np.cos(n_t)])
corner2 = -n_arrow_direction * (radius_max + arrow_length)
name = "arc" + str(i)
arc_mesh = self._create_arc(
pointa=corner1,
pointb=corner2,
center=(0, 0, 0),
resolution=100,
negative=False,
name=name,
color=line_color,
)
self.all_scene_actors["annotations"]["incident_wave"][name] = arc_mesh
return True
@pyaedt_function_handler()
@graphics_required
def _create_arrow(self, start, direction, scale, name, color):
arrow = pv.Arrow(start=start, direction=direction, scale=scale)
arrow_object = SceneMeshObject()
arrow_object.name = name
arrow_object.color = color
arrow_object.line_width = 5
arrow_object.specular = 0.5
arrow_object.lighting = True
arrow_object.smooth_shading = False
arrow_object.mesh = arrow
return MeshObjectPlot(arrow_object, arrow_object.get_mesh())
@pyaedt_function_handler()
@graphics_required
def _create_arc(self, pointa, pointb, center, resolution, negative, name, color):
arc = pv.CircularArc(pointa=pointa, pointb=pointb, center=center, resolution=resolution, negative=negative)
arc_object = SceneMeshObject()
arc_object.name = name
arc_object.color = color
arc_object.line_width = 5
arc_object.specular = 0.5
arc_object.lighting = True
arc_object.smooth_shading = False
arc_object.mesh = arc
return MeshObjectPlot(arc_object, arc_object.get_mesh())
@pyaedt_function_handler()
@graphics_required
def _create_cone(self, center, direction, radius, height, resolution, name, color):
cone = pv.Cone(center=center, direction=direction, radius=radius, height=height, resolution=resolution)
cone_object = SceneMeshObject()
cone_object.name = name
cone_object.color = color
cone_object.line_width = 5
cone_object.specular = 0.5
cone_object.lighting = True
cone_object.smooth_shading = False
cone_object.mesh = cone
return MeshObjectPlot(cone_object, cone_object.get_mesh())
@pyaedt_function_handler()
@graphics_required
def _create_line(self, pointa, pointb, name, color):
line = pv.Line(pointa=pointa, pointb=pointb)
line_object = SceneMeshObject()
line_object.name = name
line_object.color = color
line_object.line_width = 5
line_object.mesh = line
return MeshObjectPlot(line_object, line_object.get_mesh())
@pyaedt_function_handler()
@graphics_required
def __get_pyvista_range_profile_actor(
self,
xpos,
ypos,
plot_data,
cpos,
axis_vec,
plot_type="Line",
data_conversion_function="",
scene_actors=None,
extents=None,
):
"""Get PyVista actor for range profile 3D scene."""
if extents is None: # pragma: no cover
extents = [0, 10, 0, 10, 0, 10]
plot_type_lower = plot_type.casefold()
actor = None
if not plot_type_lower == "projection":
xyz_pos = np.stack((xpos, ypos, plot_data)).T
actor = pv.lines_from_points(xyz_pos)
actor[data_conversion_function] = cpos
if plot_type_lower in ["ribbon", "plane v", "plane h"]:
extents_radius = max(
(extents[0] ** 2 + extents[2] ** 2 + extents[4] ** 2) ** 0.5,
(extents[1] ** 2 + extents[3] ** 2 + extents[5] ** 2) ** 0.5,
)
geo_width = 2 * extents_radius
if plot_type_lower == "plane v":
norm_vect = axis_vec[1]
else:
norm_vect = axis_vec[2]
actor = actor.ribbon(width=geo_width / 2, normal=norm_vect)
elif plot_type_lower == "rotated":
v = axis_vec[0]
v_hat = v / np.linalg.norm(v)
actor.extrude_rotate(rotation_axis=v_hat, capping=True, inplace=True)
elif plot_type_lower == "extruded":
radius_max = max(
(extents[0] ** 2 + extents[2] ** 2 + extents[4] ** 2) ** 0.5,
(extents[1] ** 2 + extents[3] ** 2 + extents[5] ** 2) ** 0.5,
)
width_max = (
(actor.bounds[1] - actor.bounds[0]) ** 2
+ (actor.bounds[3] - actor.bounds[2]) ** 2
+ (actor.bounds[5] - actor.bounds[4]) ** 2
) ** 0.5
center = np.array(actor.center)
center = center - axis_vec[2] * radius_max
plane = pv.Plane(
center=center,
direction=-axis_vec[2],
i_size=width_max,
j_size=width_max,
)
actor.extrude_trim(-axis_vec[2], plane, inplace=True)
elif plot_type_lower == "projection":
if scene_actors is None: # pragma: no cover
return None
actor = pv.PolyData()
for model_actor in scene_actors.values():
mesh = model_actor.custom_object.get_mesh()
xypoints = mesh.points
xpos_ypos = np.column_stack((xpos, ypos, plot_data))
all_indices = self.__find_nearest_neighbors(xpos_ypos, xypoints)
mag_for_color = np.ndarray.flatten(cpos[all_indices])
if not mesh.__class__.__name__ == "PolyData": # pragma: no cover
mesh_triangulated = mesh.triangulate()
model_actor.custom_object.mesh = pv.PolyData(mesh_triangulated.points, mesh_triangulated.cells)
else:
model_actor.custom_object.mesh.clear_data()
model_actor.custom_object.mesh[data_conversion_function] = mag_for_color
actor += model_actor.custom_object.mesh
else: # pragma: no cover
raise ValueError(f"Invalid plot type: {plot_type}.")
return actor
[docs]
@staticmethod
def stretch_data(data, scaling_factor, offset):
"""
Stretches and scales the input data to a specified range.
This method normalizes the input data between its minimum and maximum values and then applies
a linear transformation using the formula: ``scaled_data = (data - min) / (max - min) * m + b``.
The parameters ``m`` and ``b`` control the scaling and shifting of the normalized data.
Parameters
----------
data : numpy.ndarray or pandas.Series
The input data array or series to be stretched.
scaling_factor : float
The scaling factor applied to the normalized data.
offset : float
The offset added to the scaled data after normalization.
Returns
-------
numpy.ndarray or pandas.Series
Transformed data
Examples
--------
>>> data = np.array([1, 2, 3, 4, 5])
>>> stretched_data = stretch_data(data, 2, 1)
>>> print(stretched_data)
[1. 1.5 2. 2.5 3. ]
"""
return (data - data.min()) / (data.max() - data.min()) * scaling_factor + offset
@staticmethod
def __find_nearest_neighbors(orig_points, target_points, grid_origin=None, u=None, v=None, w=None, grid_shape=None):
"""
Fast nearest neighbor search for points on a 2D or 3D grid in 3D, using grid geometry.
If grid_origin, u, v, (optionally w), and grid_shape are provided, uses projection-based search.
Otherwise, falls back to brute-force search.
Parameters
----------
orig_points : numpy.ndarray
Array of shape (n, 3) representing the original grid points (x, y, z).
target_points : numpy.ndarray
Array of shape (m, 3) representing the target points (x, y, z).
grid_origin : numpy.ndarray, default: ``None``
3D origin of the grid (shape (3,)).
u : numpy.ndarray, default: ``None``
3D basis vector for grid axis 0 (shape (3,)).
v : numpy.ndarray, default: ``None``
3D basis vector for grid axis 1 (shape (3,)).
w : numpy.ndarray, default: ``None``
3D basis vector for grid axis 2 (shape (3,)), for 3D grids.
grid_shape : tuple, default: ``None``
(n_i, n_j) or (n_i, n_j, n_k) shape of the grid.
Returns
-------
numpy.ndarray
Array of shape (m,) containing the indices of the
nearest neighbors in the ``orig_points`` parameter
for each target point.
"""
if grid_origin is not None and u is not None and v is not None and grid_shape is not None:
if w is not None and len(grid_shape) == 3:
# 3D grid
u_ = np.column_stack((u, v, w)) # 3x3
u_pinv = np.linalg.pinv(u_) # 3x3
n_i, n_j, n_k = grid_shape
indices = []
for p in target_points:
ijk = u_pinv @ (p - grid_origin)
i, j, k = np.round(ijk).astype(int)
i = i % n_i
j = j % n_j
k = k % n_k
idx = i * (n_j * n_k) + j * n_k + k # flatten index
indices.append(idx)
return np.array(indices)
else:
# 2D grid (original logic)
n = np.cross(u, v)
n = n / np.linalg.norm(n)
u_ = np.column_stack((u, v)) # 3x2
u_pinv = np.linalg.pinv(u_) # 2x3
n_i, n_j = grid_shape
indices = []
for p in target_points:
# Project onto grid plane
p_proj = p - np.dot(p - grid_origin, n) * n
# Get (i, j) coordinates
ij = u_pinv @ (p_proj - grid_origin)
i, j = np.round(ij).astype(int)
i = i % n_i
j = j % n_j
idx = i * n_j + j # flatten index
indices.append(idx)
return np.array(indices)
else: # pragma: no cover
# Fallback: brute-force search
distances = ((orig_points[:, np.newaxis] - target_points) ** 2).sum(axis=2)
return np.argmin(distances, axis=0)
@staticmethod
def __add_mesh(mesh_object, plotter, mesh_type="results"):
"""Add a mesh to the plotter with additional options."""
options = {}
if getattr(mesh_object, "custom_object", None):
if mesh_type == "model":
options = mesh_object.custom_object.get_model_options()
elif mesh_type == "annotations":
options = mesh_object.custom_object.get_annotation_options()
else:
options = mesh_object.custom_object.get_result_options()
plotter.plot(mesh_object.custom_object.get_mesh(), **options)
@pyaedt_function_handler()
def __get_model_extent(self):
"""
Calculate the 3D extent of the model by evaluating the bounding box dimensions of each mesh object in the scene.
This method retrieves the maximum and minimum coordinates in the x, y, and z
directions for all mesh objects stored under the "model" key in `self.all_scene_actors`.
The bounding box of each mesh is assessed, and the overall bounds for the entire
model are determined by taking the min/max values from these individual bounding boxes.
"""
x_max, x_min, y_max, y_min, z_max, z_min = [], [], [], [], [], []
if len(self.all_scene_actors["model"]) == 0:
x_max = [1]
x_min = [-1]
y_max = [1]
y_min = [-1]
z_max = [1]
z_min = [-1]
for each in self.all_scene_actors["model"].values():
b = each.custom_object.get_mesh().bounds
x_max.append(b[1])
x_min.append(b[0])
y_max.append(b[3])
y_min.append(b[2])
z_max.append(b[5])
z_min.append(b[4])
self.__x_max, self.__x_min = max(x_max), min(x_min)
self.__y_max, self.__y_min = max(y_max), min(y_min)
self.__z_max, self.__z_min = max(z_max), min(z_min)
@pyaedt_function_handler()
@graphics_required
def __get_geometry(self):
"""Get 3D meshes."""
model_info = self.model_info
obj_meshes = {}
first_value = next(iter(model_info.values()))
self.__model_units = first_value[3]
for object_in in model_info.values():
relative_cad_path, color, opacity, units = object_in
relative_path = Path(relative_cad_path)
name = relative_path.stem
relative_path = Path("geometry") / relative_path
cad_path = Path(self.rcs_data.output_dir) / relative_path
try:
conv = AEDT_UNITS["Length"][units]
except Exception: # pragma: no cover
conv = 1
if cad_path.exists():
mesh = pv.read(str(cad_path))
mesh.scale(conv)
else: # pragma: no cover
self.__logger.warning(f"{cad_path} does not exist.")
return False
color_cad = [i / 255 for i in color]
model_object = SceneMeshObject()
model_object.color = color_cad
model_object.opacity = opacity
model_object.name = name
model_object.mesh = mesh
mesh_object = MeshObjectPlot(model_object, model_object.get_mesh())
obj_meshes[model_object.name] = mesh_object
return obj_meshes
[docs]
@staticmethod
def rotate_point(x1, y1, z1, rotation):
"""
Rotate a set of 3D points by the specified azimuth, elevation, and twist angles.
Parameters
----------
x1 : np.ndarray
Array of x-coordinates of the points to be rotated.
y1 : np.ndarray
Array of y-coordinates of the points to be rotated.
z1 : np.ndarray
Array of z-coordinates of the points to be rotated.
rotation : dict
Dictionary containing rotation angles in degrees with the following keys:
- "azimuth": Rotation angle around the z-axis.
- "elevation": Rotation angle around the y-axis.
- "twist": Rotation angle around the x-axis.
Returns
-------
x2 : np.ndarray
Array of rotated x-coordinates.
y2 : np.ndarray
Array of rotated y-coordinates.
z2 : np.ndarray
Array of rotated z-coordinates.
Notes
-----
The rotation is performed in the following order:
1. Azimuth (z-axis)
2. Elevation (y-axis)
3. Twist (x-axis)
"""
az_rad = np.radians(rotation["azimuth"])
el_rad = np.radians(rotation["elevation"])
tw_rad = np.radians(rotation["twist"])
# Rotation around z-axis (azimuth direction)
maz = np.array([[np.cos(az_rad), -np.sin(az_rad), 0], [np.sin(az_rad), np.cos(az_rad), 0], [0, 0, 1]])
# Rotation around y-axis (elevation direction)
mel = np.array([[np.cos(el_rad), 0, np.sin(el_rad)], [0, 1, 0], [-np.sin(el_rad), 0, np.cos(el_rad)]])
# Rotation around x-axis (twist direction)
mtw = np.array([[1, 0, 0], [0, np.cos(tw_rad), -np.sin(tw_rad)], [0, np.sin(tw_rad), np.cos(tw_rad)]])
# Combine the rotations
r = np.dot(mtw, np.dot(mel, maz))
# Twist is applied to rotated local coordinate. This corrects to
# rotating along global coordinate.
r = uvecs2rot(r[0, :], r[1, :])
# apply the rotations
points = np.column_stack((np.ravel(x1), np.ravel(y1), np.ravel(z1)))
rotated_points = np.dot(r, points.T).T
x2 = rotated_points[:, 0].reshape(x1.shape)
y2 = rotated_points[:, 1].reshape(y1.shape)
z2 = rotated_points[:, 2].reshape(z1.shape)
return x2, y2, z2
def uvecs2rot(u, v):
"""
Obtain roll-pitch-yaw angles from two unit vectors.
It obtains the angles of rotation in term of its RPY (roll, pitch, yaw)
from the two unit vectors introduced as parameter along with their
selector
Roll-Pitch-Yaw convention follows Savant convention for defining
coordinate systems: Start with local coordinate system aligned
with global coordinate system. Then perform left-handed rotations
of the local coordinate system ABOUT ITS OWN CURRENT AXES in the
following order: Roll (xlcl, first), Pitch (ylcl, second),
Yaw (zlcl, third and last).
Args:
u (np.ndarray): First input vector. Modified in-place to be unit length.
v (np.ndarray): Second input vector. Modified in-place to be unit length.
Returns
-------
np.ndarray: A 3x3 rotation matrix with `u`, `v`, and `w` as columns.
Raises
------
ValueError: If `u` or `v` is a zero vector.
"""
u /= np.linalg.norm(u)
v /= np.linalg.norm(v)
w = np.cross(u, v)
return np.column_stack((u, v, w))
class SceneMeshObjectType(Enum):
MESH = 1
VOLUME = 2
class SceneMeshObject:
"""
A class representing a custom 3D mesh object with visualization properties.
This class defines a 3D mesh object with customizable properties.
It provides methods to retrieve the mesh, its associated rendering options, and annotation properties for
visualization in PyVista.
"""
@graphics_required
def __init__(self):
# Public
self.name = "CustomObject"
self.opacity = 1.0
self.color = None
self.color_map = "jet"
self.line_width = 1.0
self.specular = 0.0
self.lighting = False
self.smooth_shading = True
self.edge_color = None
self.show_edges = True
self.scalar_dict = dict(color="#000000", title="Dummy")
self.show = True
self.object_type = SceneMeshObjectType.MESH
self.clim = None
# Private
self.__mesh = pv.Cube()
self.__original_points = self.mesh.points.copy()
self.__z_offset = 0.0
self.__scale_factor = 1.0
self.__default_clim = None
self.__plot_type = None
self.__default_mesh = self.mesh
@property
def mesh(self):
"""Get the mesh object."""
return self.__mesh
@mesh.setter
def mesh(self, val):
self.__mesh = val
self.__original_points = val.points.copy()
if self.__default_mesh is None: # pragma: no cover
self.__default_mesh = val
@property
def z_offset(self):
"""Offset in the Z direction."""
return self.__z_offset
@z_offset.setter
def z_offset(self, val):
translation_distance = val
# Calculate the new points by applying the translation to the original points
new_points = self.__original_points.copy()
new_points[:, 2] += translation_distance # Apply Z translation
# Update the mesh with the new points
self.__mesh.points = new_points
self.__z_offset = val
@property
def scale_factor(self):
"""Get the current scale factor."""
return self.__scale_factor
@scale_factor.setter
def scale_factor(self, val):
"""Set a new scale factor and update the mesh accordingly."""
# Calculate the center of the mesh for scaling
center = self.__mesh.points.mean(axis=0) # Center of the original mesh
# Calculate the new points by scaling relative to the original points and the mesh center
new_points = center + (self.__original_points - center) * float(val)
self.__mesh.points = new_points
# Update the scale factor
self.__scale_factor = val
@property
def default_clim(self):
"""Default minimum and maximum field values."""
return self.__default_clim
@default_clim.setter
def default_clim(self, value):
self.__default_clim = value
@property
def default_mesh(self):
"""Original mesh."""
return self.__default_mesh
@default_mesh.setter
def default_mesh(self, value):
self.__default_mesh = value
@property
def plot_type(self):
"""Plot type."""
return self.__plot_type
@plot_type.setter
def plot_type(self, value):
self.__plot_type = value
def reset_scene(self):
"""Reset the mesh to its original position and size."""
self.mesh.points = self.__original_points.copy() # Restore the original points
self.__z_offset = 0.0 # Reset the Z-offset
self.__scale_factor = 1.0 # Reset the scale factor
def get_mesh(self):
"""Retrieve the mesh object.
Returns
-------
pyvista.PolyData or pyvista.UnstructuredGrid
Mesh object representing the 3D geometry.
"""
return self.mesh
def get_model_options(self):
"""Retrieve the visualization options for the mesh.
Returns
-------
dict
A dictionary with the color and opacity settings for rendering the model.
"""
return {"color": self.color, "opacity": self.opacity}
def get_annotation_options(self):
"""Retrieve the annotation options for the mesh.
Returns
-------
dict
A dictionary with the color and line width settings for annotating the model.
"""
return {
"color": self.color,
"line_width": self.line_width,
"specular": self.specular,
"lighting": self.lighting,
"smooth_shading": self.smooth_shading,
"opacity": self.opacity,
"show_edges": self.show_edges,
"edge_color": self.edge_color,
}
def get_result_options(self):
"""Retrieve the result options for the mesh.
Returns
-------
dict
A dictionary with the settings for results the model.
"""
options = {
"opacity": self.opacity,
"cmap": self.color_map,
"scalar_bar_args": self.scalar_dict,
"clim": self.clim,
}
if self.object_type != SceneMeshObjectType.VOLUME:
options["line_width"] = self.line_width
if self.color:
options["color"] = self.color
return options
