Plotting Dynamic Contours on a Floor plan

Hello,

So i have this function and want to make the contour lines dynamic/adaptable, meaning that i want the contours to recognise that there is a wall/boundary and from there to use a different equation. How would i do that? it seems i hit a brick wall and can not think of the solution.

Thank you in advance!

NOTE: The image i load is a floor plan.

% Function to detect walls in the image and plot them in 2D
function visualizeWalls2D(lengthField, widthField, distanceEditField, sourceDropdown, sourceData, doseEditField, rateEditField, activityEditField, density, resultTable)
%Retrieve source data
selectedSource = sourceData.(sourceDropdown.Value);
materials = fieldnames(selectedSource.TVLe);

% Retrieve the image from the application data
img = getappdata(0, 'FloorPlanImage');
if isempty(img)
    return;
end

% Get the user-defined Length, Width, and Thickness
floorLength = lengthField.Value;
floorWidth = widthField.Value;

if isempty(floorLength) || isempty(floorWidth)
    uialert(gcf, 'Please enter valid dimensions!', 'Invalid Input');
    return;
end

% Convert image to grayscale and adjust contrast
gs = im2gray(img);
gs = imadjust(gs);

% Detect edges using edge detection (e.g., Canny method)
edges = edge(gs, 'Canny', [0.1 0.3]);

% Find the boundaries of the objects in the binary image
[B, ~] = bwboundaries(edges, 'noholes');

for i = 1:length(materials)
    material = materials{i};
    % Create a new figure for the 2D visualization
    fig = figure("Name", "2D " + material + " Wall Visualization");
    ax = axes('Parent', fig); % Create axes in the new figure
    view(ax, 2); % Set view to 2D
    xlabel(ax, 'X');
    ylabel(ax, 'Y');

    hold(ax, 'on'); % Hold the plot for multiple walls
    for k = 1:length(B)
        boundary = B{k};

        % Normalize the boundary coordinates to be between 0 and 1
        normalized_x = boundary(:,2) / size(edges, 2); % Normalize by image width
        normalized_y = 1 - (boundary(:,1) / size(edges, 1)); % Normalize by image height

        % Scale normalized coordinates based on the user input
        scaled_x = normalized_x * floorLength;
        scaled_y = normalized_y * floorWidth;

        % Plot the wall with user-defined thickness
        plot(ax, scaled_x, scaled_y, 'k', 'LineWidth', 2); % Bottom boundaries
        hold on
    end
    distance = zeros(1,4);
    sourceX=distanceEditField{1}.Value;
    sourceY=distanceEditField{2}.Value;
    plot(ax, sourceX, sourceY, 'r*');
    hold on
    for j = 1:4
        distance(j) = distanceEditField{j}.Value;
    end
    [maxdistance, idx] = max(distance);
    [d, theta] = meshgrid(linspace(0.25, maxdistance, 200), linspace(0, 2*pi, 200));
    dosedist = zeros(size(d));
    for l = 1:numel(d)
        [row, col] = ind2sub(size(d),l);
        distFromSource = d(row, col);
        xPos = sourceX + distFromSource * cos(theta(row, col));
        yPos = sourceY + distFromSource * sin(theta(row, col));
        distanceToSource = sqrt((xPos - scaled_x).^2 + (yPos - scaled_y).^2);
        if distFromSource >= abs(distanceToSource)
            % Apply attenuation if within wall boundaries
            dosedist(row, col) = (activityEditField.Value * selectedSource.RAKR * doseEditField.Value) ./ (distFromSource.^2 * rateEditField.Value * 60) * exp(-selectedSource.(material) * density.(material) * resultTable.Data{j, 2*idx});
        else
            % Normal dose calculation without attenuation
            dosedist(row, col) = (activityEditField.Value * selectedSource.RAKR * doseEditField.Value) ./ (distFromSource.^2 * rateEditField.Value * 60);
        end
    end
    maxDose = max(dosedist(:), [], 'omitnan');
    contourLevels = logspace(log10(maxDose/100), log10(maxDose), 20); % 20 contour levels with closer spacing near lower doses
    contour(sourceX + d.*cos(theta), sourceY + d.*sin(theta), dosedist, contourLevels, 'ShowText', 'on');
    colorbar;
    hold (ax, 'off')
end
end

[Aquatris]

Can you provide a runable code and show which parts in the contour you want things to be different?

[ALEXANDROS]

This is an image of an example from the code. In the programme i load the floor plan image, detect the boundaries with Canny and then i calculate the dose values for the maximum distance out of the four that the user had defined, plotting them on the floor plan as contour lines. But it seems that the programme does not take into consideration the boundaries when plotting the contours. It should exponentially decrease at the wall, thus resulting in non-uniform/non-concentric distribution.

I will provide a runable code later on as there is a lot of code and much of it is for the GUI i created and i am also fairly new to matlab, it might take me a while. What you see is a month or two of a lot of trial and error.

Answer 1

[Shishir Reddy]

Hi Alexandros,

From the provided code, I understand that you have the boundaries of the walls of a floor plan, extracted using MATLAB's bwboundaries function. The dose calculation is based on user-defined parameters such as source position, activity, and material properties.

As per my understanding, you want to implement a method to check if a line from the radiation source to any point in the dose distribution grid intersects with any of the detected wall boundaries. Additionally, you want to modify the dose calculation logic to account for these wall intersections. To achieve this in MATLAB, please consider the following steps:

1. Intersection Check Function: You need to create a function to check if a line intersects with any boundary. This function should take source coordinates, target coordinates, and the boundary data as input, and return whether an intersection occurs.

function intersects = checkIntersection(source, target, boundaries)
    intersects = false;
    for k = 1:length(boundaries)
        boundary = boundaries{k};
        for n = 1:size(boundary, 1) - 1
            if isIntersecting(source, target, boundary(n,:), boundary(n+1,:))
                intersects = true;
                return;
            end
        end
    end
end

Here, isIntersecting is another function you need to create to check whether two lines intersect. This can be achieved using a simple cross product.

2. Dose Calculation: Apply conditional logic to different equations based on the presence or absence of walls.

if checkIntersection([sourceX, sourceY], [xPos, yPos], B)
    % Attenuation
    dosedist(row, col) = baseDose * exp(-attenuationFactor);
else
    dosedist(row, col) = baseDose;
end

These changes introduce a mechanism to check for intersections with walls and adjust the dose calculation accordingly. Please note that these are just snippets, so ensure you integrate them into your original code structure.

I hope this helps.

[ALEXANDROS]

Hello, sorry for the late reply I have been occupied with another section of my code. I will check it out as soon as I can, thank you very much for your reply.

Answer 2

[Michelle]

Hi Alexandos,

If you want the contour line to remain circular with their color changing/attenuating after meeting a wall, contour plot is not adapted to what you want to do. Indeed, the line will follow the dose value, meaning it will break it circular shape when meeting a wall such like this:

I recommend plotting the contour with non-attenuated dose values in order to keep the nice circular waves and using a white surf plot with dose-modulates transparency (AlphaData) on top of the contour lines to attenuate their color:

% plot contour lines using non-attenuated dose for line color/level, dosedist_0
contour(sourceX + d.*cos(theta), sourceY + d.*sin(theta), dosedist_0/maxDose, contourLevels);

% plot white surface using attenuated dose for transparency, dosedist
surface(sourceX + d.*cos(theta), sourceY + d.*sin(theta), zeros(size(theta)), 'alphadata',...
    1-dosedist/maxDose, 'edgecolor', 'none', 'facealpha', 'flat', 'facecolor', [1,1,1]);

% then plot your labyrinth on top of all
for k = 1:length(B)
    
    % [...]
    
    % Plot the wall with user-defined thickness
    plot(ax, scaled_x, scaled_y, 'k', 'LineWidth', 2); % Bottom boundaries
    hold on
end

It will give something like this:

Or with uniformly black contour levels:

% plot contour lines using non-attenuated dose for line color/level, dosedist_0
contour(sourceX + d.*cos(theta), sourceY + d.*sin(theta), dosedist_0/maxDose, contourLevels, ...
    'edgecolor', [0,0,0]);

EDIT:

I understand from your comment that your question is concerning rather the dose calculation than the plot.

Here is the code I used:

% Convert boundaries coordinates (scaled_xy) to polar (rwall,twall)
xwall = scaled_x - sourceX;
ywall = scaled_y - sourceY;
rwall = sqrt(xwall.^2 + ywall.^2);
twall = atan(ywall./xwall);
twall(xwall>0 & ywall<0) = twall(xwall>0 & ywall<0)+2*pi;
twall(xwall<0) = twall(xwall<0)+pi;
twall(xwall==0 & ywall>0) = pi/2;
twall(xwall==0 & ywall<0) = 3*pi/2;

% Sample space around source
Nt = 200;
Nr = 200;
t = linspace(0, 2*pi, Nt);
r = linspace(0.25, maxdistance, Nr);
[d, theta] = meshgrid(r, t);

% Normal dose calculation without attenuation
dosedist_0 = repmat((activityEditField.Value * selectedSource.RAKR * ...
    doseEditField.Value) ./ (r.^2 * rateEditField.Value * 60),[Nt,1]);
    
% Attenuation factor
fact = exp(-selectedSource.(material) * density.(material) * ...
    resultTable.Data{j, 2*idx});

% Use sine and cosine
coswall = cos(twall);
sinwall = sin(twall);
cost = cos(t);
sint = sin(t);

dosedist = dosedist_0;
for n = 1:Nt
    
    % Find all places where the dose meets boundaries 
    nwall = find(...
            (...
                ((coswall(1:end-1)>=cost(n) & coswall(2:end)<=cost(n)) | ...
                 (coswall(2:end)>=cost(n) & coswall(1:end-1)<=cost(n)) ...
                ) ...
                & sign(sinwall(1:end-1))==sign(sint(n)) ...
                & sign(sinwall(2:end))==sign(sint(n))...
            ) | (...
                ((sinwall(1:end-1)>=sint(n) & sinwall(2:end)<=sint(n)) | ...
                 (sinwall(2:end)>=sint(n) & sinwall(1:end-1)<=sint(n)) ...
                ) ...
                & sign(coswall(1:end-1))==sign(cost(n)) ...
                & sign(coswall(2:end))==sign(cost(n))...
            ) ...
        );
    
    % Attenuate dose every time it meets a new boundary
    [~,ord] = sort(rwall(nwall));
    nwall = nwall(ord);
    for m = nwall'
        dosedist(n,r>=mean(rwall([m,m+1]))) = dosedist(n,r>=mean(rwall([m,m+1])))*fact;
    end
end

% Plot attenuated dose with contours
maxDose = max(dosedist_0(:), [], 'omitnan');
contourLevels = logspace(-2, 0, 20); % 20 contour levels with closer spacing near lower doses
contour(sourceX + d.*cos(theta), sourceY + d.*sin(theta), dosedist/maxDose, ...
    contourLevels);

I hope it works for you.

[ALEXANDROS]

Hello, sorry for the late reply I have been occupied with another section of my code. The first image is actually what I really want. What code did you use to accomplish that? Thank you very much for your reply.

[Michelle]

Hi Alexandros, see the edit in my answer. Please tell me if it works for you.