Ray Tracing with raytrace(): define a reception local area instead of a single point

Hi,

I'm trying to use Matlab's raytrace() function from the Communications toolbox with the Shooting & Bouncing Rays method to compute the received rays in a local area around the receiver RX, instead of just in a single point (RX).

The context is the following: I'm trying to compute reflections of RF waves on a human body, using Matlab's raytrace() function. I model the human body as a cylinder, compute a triangulation (needed for the raytrace() function) for it and feed that to the raytrace() function. Let's say I put the cylinder at the perfect angle w.r.t. the transmitter TX and the receiver RX, a ray corresponding to a single reflection is found. However if I move the cylinder of 10 cm the reflection is not found, which is logical. It is illustrated on the plot below (left: a reflection is found, right: no reflection is found):

However since the human body is not a perfect cylinder, in reality there is a high chance that for such a small displacement (10 cm) from the ideal position a reflection would still be present. That is why I would like to compute all rays arriving in a local area around RX as it is usually done in ray tracing simulations to compute the received power at RX, rather than at the RX precise point. This way I would detect rays even when the human body is not at the perfect position w.r.t. TX and RX.

I've already considered instantiating several RXs in a local area, but then I detect several rays while there should be only one, so if I add a new cylinder in the scene then I will not be able to discriminate those artificial rays from the ones coming from that new cylinder. I've also dug in the raytrace() function code. At line 405 it calls another raytrace() function, which at line 313-323 makes the test I want to perform. It defines a radius around RX and checks if the ray falls within that radius (I'm using Matlab 2021a):

% Calculate reception radius, i.e., alpha*d/sqrt(3) 
radius = ((rayDist + srcPrjDist) .* maxAngularSep)/sqrt(3);

% Detect which rays hit the Rx and mitigate double counting.
% Because the first segment for LOS is already covered, no need to
% perform detection for it.
if refIdx > 1
    % Flags of the rays that reach Rx. Note that intDist can be NaN
    % if the ray hits a sharp edge. But (NaN > a) = 0 for any value
    % of a.
    isRxReached = (srcPrjDist > 0) & (rxPrjDist <= radius) & ... 
        (intDist >= srcPrjDist);  

A way to proceed would be to extend that radius, but the function cannot be modified unfortunately, and shadowing the function with a modified local copy is usually bad practice (and doesn't work in this case since some validations are done within the function, that require code accessible only from the toolbox folder).

So does someone have an idea on how to proceed to compute rays in a local area? I include herebelow the code I use to generate TX, RX, the cylinder and the propagation model.

%% TX/RX
fc = 2.45e9;
tx = txsite("cartesian",  "AntennaPosition", [0; 0; 1],"TransmitterFrequency", fc);
rx = rxsite("cartesian","AntennaPosition", [0; 3; 1]);
%figure;
scatter3(tx.AntennaPosition(1,:), tx.AntennaPosition(2,:), tx.AntennaPosition(3,:), 'sr', 'filled');hold on;
scatter3(rx.AntennaPosition(1,:), rx.AntennaPosition(2,:), rx.AntennaPosition(3,:), 'sb', 'filled');hold on;
xlabel('x');ylabel('y');zlabel('z');
legend('TX','RX');

%% Model human body as a cylinder
r = 1; % [m]
height = 1.75; % [m]
theta = linspace(0,2*pi,36).';
xc = 1.5; % center of the cylinder
yc = 1.6;
x = repmat(xc + r*cos(theta),[2 1]);
y = repmat(yc + r*sin(theta),[2 1]);
z = [0*ones(length(theta),1);height*ones(length(theta),1)];
P = [x y z];
[k,vol] = convhulln(P); % convex hull: set of points defining the boundary of the volume
trisurf(k,P(:,1),P(:,2),P(:,3),'FaceColor','magenta')
TR = triangulation(k,P);

%% Propagation model and raytracing
pm = propagationModel("raytracing","CoordinateSystem","cartesian","Method","sbr","MaxNumReflections",2,"SurfaceMaterial","perfect-reflector",...
    "AngularSeparation","low");
rays = raytrace(tx, rx, pm, 'Map', TR); 
noMPC = 1;
rays_n = rays{1};
if length(rays_n) > 1
    refl = rays_n(2).ReflectionLocations; % the first element of rays_n is the direct path
    line1 = [tx.AntennaPosition refl];
    line2 = [refl rx.AntennaPosition];
    plot3(line1(1,:),line1(2,:),line1(3,:),'b');
    plot3(line2(1,:),line2(2,:),line2(3,:),'b');
    noMPC = 0;
else
    disp('No MPC, only direct path'); 
end

Thank you for your help!

Answer 1

[Yue Shang]

Hi Laurent,

First, creating a finer mesh for the cylinder may help. For example, you will find one ray in the failing case if you define

theta = linspace(0,2*pi,360).';

I would suggest instantiating multiple rx sites in the local area as you mentioned, instead of modifying the reception sphere radius which can invalidate the shooting bouncing rays (SBR) method. For example, you can define 27 rx sites in a cubic around the primary rx site of interest:

[x, y,z] = meshgrid(-0.1:0.1:0.1, 2.9:0.1:3.1, 0.9:0.1:1.1);
rx = rxsite("cartesian","AntennaPosition", [x(:),y(:),z(:)]');

In the ray tracing result, rays{14} would be the "true" rays arriving at the primary rx and the other cell array elements contain the "artificial" rays arriving at the local area, for any number of cylinders in the scene.

[Laurent Storrer]

Hi Yue,

Thanks a lot for your answer! Indeed the finer mesh helps a bit but only for a small displacement w.r.t. the ideal position, so I'm not sure that it is worth the extra computational burden (I plan to have a lot of cylinders, around 50).

What do you mean by "invalidate the shooting bouncing rays" if the radius is increased? If this radius is defined in meters, increasing it (by a small factor, e.g. 1.5) would just be equivalent to considering a larger reception area, wouldn't it?

I tried instantiating 27 RX sites in a cube as you suggested, and indeed it allows to have rays for a wider range of positions of the cylinder. However I'm not sure it fulfills my needs, because what I would like to obtain is one unique ray but in an increased reception area: an "artificial" ray would still carry a valuable information event if rays{14} (the "true" ray) doesn't exist, because the artificial ray would correspond to a part of the human body that is not modelled by the cylinder and that can still create a reflection. But with several RXs I would have several of those rays, while I would like to have only one ray, "true" or "artificial", corresponding to the target. (For information, this is for a bistatic radar application, in which a target will often give a reflection even if it is not located at the perfect position w.r.t. TX and RX because it has a body shape quite different from a cylinder, with different parts, orientations etc.)

• So wouldn't increasing the reception radius be equivalent to instantiating several RX, without suffering from the problem of having multiple rays for one cylinder and with less computation time since there is only one RX?

Thank you again for your help.

Best regards,

Laurent

[Yue Shang]

Hi Laurent,

In the SBR ray tracing method, each ray represents a cone-shape tube whose size grows as the ray travels. The overlap between adjacent tubes should be minimal, for which the current reception sphere radius is designed. Increasing the reception sphere will likely end up with more than one ray for one rx, while there should be only one. You can refer to the following page for some description of this ray tracing method:

https://www.mathworks.com/help/comm/ug/choose-a-propagation-model.html

When you have multiple "artificial" rays, you may just pick one from them for your applications.

Hope this helps.

[Laurent Storrer]

Hi Yue,

Ok I see so increasing the radius would give rise to false rays.

I will go for your suggested approach with the multiple RXs in a cube then, and I think I will also try other surfaces than a cylinder to model the human body.

Thank you for your help!

Best regards,

Laurent

Answer 2

[Saad Humayun]

Would you by any chance know how to access this part of raytrace(0 function in the newer versions of MATLAB? I know this is super late, but I would like to see the calculation of reception sphere and flagging logic similarly in the later versions too, as some of my rays in simulations are being missed and I want to be able to explain that phenomenon.

[Laurent Storrer]

Hi Saad,

If you copy the code in my original question, you can right-click on the raytrace function and select "open raytrace". Then, inside that function, there is a call to comm.internal.raytrace, at line 563 in Matlab R2023b. It is a private function that can only be called from the first raytrace function. To be able to open it, you need to put a breakpoint at the corresponding line in the first raytrace function (so at line 563). In that function comm.internal.raytrace you can check the details of the ray computation. But from what I've just seen, it seems to call a C++ function for the actual ray computation. To find this one, maybe you will need to check manually inside Matlab's directories.

Hope this helps,

Laurent

[Saad Humayun]

I was just exploring comm.internal.raytrace in R2023b, as that is the one I was using earlier. From what you said, that is literally the best conclusion I came up with, because unfortunately I wasn't lucky enough to find any actual flagging or computation, just that C++ line instead which obviously I don't think users have access to.

Thanks a lot Laurent, super grateful for the comment!

[Saad Humayun]

Aside from this, would you think logically and theoretically MATLAB would keep the calculations similar even for the later versions when it comes to reception mechanisms of rays and receiver sites? Even though the exact lines of code aren't visible in these versions.