Issues plotting R,theta, T from X,Y,Z data

Hello everyone!

I have data of a cylinder with certain tempretures at each point and the data is unfortunatley in X,Y,Z coordinates (as if the cylinder is floating in a cartesian system). I wrote code that transforms th data to R, theta coordinates as well as pick the outer and inner radii and plot them on a 2d plane (theta vs X with T colour map). The problem I have is that the data from around 0-0.6 radians is messed up with either missing or wrong points plotted. I do use a tolerance for matlab to check if the point being checked is outer or inner but increasing the tolerance anymore makes the graph intefer with inner layers. Here is the code below:

clc;
clear;

coords = readtable("Al_Case1.xlsx");
coords = table2array(coords);
X = coords(:,1);
Y = abs(coords(:,2));
Z = coords(:,3);
T = coords(:,4);
Z_center = mean(Z);
Y_center = mean(Y);
R = ((Z-Z_center).^2 + (Y-Y_center).^2).^0.5;
Theta = atan2(Z-Z_center, Y-Y_center);
itt = 0;
itti = 0;
Ro = max(R);
Ri = min(R);

for i = 1:length(R)
    if Theta(i) < 0
        Theta(i) = Theta(i)+pi;
    end
    if Ro-R(i) <= 0.25 || R(i) > Ro
        if itt >= 1 && ismember(Theta(i),Thetao) && X(i) == thicknesso(ismember(Theta(i),Thetao))
            io = ismember(Theta(i),Thetao);
            if R(i) > Router(io)
                Router(i) = R(i);
                To(i) = T(i);
                thicknesso(i) = X(i);
                Thetao(i) = Theta(i);
                itt = itt+1;
            end
        else
            Router(i) = R(i);
            To(i) = T(i);
            thicknesso(i) = X(i);
            Thetao(i) = Theta(i);
            itt = itt+1;
        end
        
    end
    if R(i)-Ri <= 0.25 || R(i) < Ri
        if itti >= 1 && ismember(Theta(i),Thetai) && X(i) ==thicknessi(ismember(Theta(i),Thetai)) 
            ii = ismember(Theta(i),Thetai);
            if R(i) < Rinner(ii)
                Rinner(i) = R(i);
                Ti(i) = T(i);
                thicknessi(i) = X(i)-300;
                Thetai(i) = Theta(i);
                itti = itti+1;
            end
        else
            Rinner(i) = R(i);
            Ti(i) = T(i);
            thicknessi(i) = X(i)-300;
            Thetai(i) = Theta(i);
            itti = itti+1;
        end
    end
end

valid_outer = Router ~= 0;
Router = Router(valid_outer);
To = To(valid_outer);
Thetao = Thetao(valid_outer);
thicknesso = thicknesso(valid_outer);

valid_inner = Rinner ~= 0;
Rinner = Rinner(valid_inner);
Ti = Ti(valid_inner);
Thetai = Thetai(valid_inner);
thicknessi = thicknessi(valid_inner);

[Thetao, idx_o] = sort(Thetao, 'descend');
thicknesso = thicknesso(idx_o);
Router = Router(idx_o);
To = To(idx_o);
Thetao = Thetao *360/pi;

[Thetai, idx_i] = sort(Thetai, 'descend');
thicknessi = thicknessi(idx_i);
Rinner = Rinner(idx_i);
Ti = Ti(idx_i);
Thetai = Thetai *360/pi;


scatter(Thetao,thicknesso,20,To,"filled")
hold on
scatter(Thetai,thicknessi,20,Ti,"filled")
colorbar;
title('Cylinder Temperature');
grid on;
hold off

I tried with some i statments to make it check if the point was already plotted and if it was to pick the point that is outermost (for plotting the outter surface) and innermost (for inner surface) but that didn't seem to change much.

Any clue or way I could go to fix it would be greatly appreciated thanks!

[Ahmed]

[Mathieu NOE]

hello

seems I have to ask every time - please share the data that goes with the code.... otherwise we're not in position to help efficiently -

[Star Strider]

Consider using the cart2sph or cart2pol function for the conversion.

[Ahmed]

Sorry for not attaching the data. Here it is.

[Ahmed]

cart2pol didn't change much

[Mathieu NOE]

sure, your code is pretty much what cart2pol function does

Answer 1

[Mathieu NOE]

hello

when plotting the X,Y,Z data you get a cylinder or something like a cylinder with elliptical shape

I removed the X dimensions and focused on ploting R vs theta and we get a mean value (Rmean = 99.2208) + a sine wave + noise or outliers

NB the temperature T is color coded (values are in the colorbar) - the dots are filled with constant size but you could also make the dots size according to T

then we can draw the upper and lower envelope (green / red) and also a smoothed version of the data (black)

question what data do you want to extract now ? remove ? I figure out maybe not correctly that the fluctuation of R is what you want , maybe after we managed to remove most of the outliers we see below the sine curve ?

code :

coords = readtable("Al_Case1.xlsx");
coords = table2array(coords);
X = coords(:,1);
Y = coords(:,2);
Z = coords(:,3);
T = coords(:,4);
Z_center = mean(Z);
Y_center = mean(Y);
R = ((Z-Z_center).^2 + (Y-Y_center).^2).^0.5;
Theta = atan2(Z-Z_center, Y-Y_center);

% keep unique values 
[Theta,ia,ib] = unique(Theta);
R = R(ia);
T = T(ia);
Rmean = mean(R);

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

% let's do some processing
Rs = smoothdata(R,'gaussian',2500);
[upperEnv,lowerEnv] = envelope(R, 100, 'peaks');

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rs,'k','linewidth',3);
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
plot(Theta,upperEnv,'r','linewidth',3);
plot(Theta,lowerEnv,'g','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

[Mathieu NOE]

with some further processing , we can remove all the dots that are under the sine wave

and do a sine fit to get a clean sine wave with T coded with color

now I assume it should be easy to use the sinewave parameters for your own needs

coords = readtable("Al_Case1.xlsx");
coords = table2array(coords);
X = coords(:,1);
Y = coords(:,2);
Z = coords(:,3);
T = coords(:,4);
Z_center = mean(Z);
Y_center = mean(Y);
R = ((Z-Z_center).^2 + (Y-Y_center).^2).^0.5;
Theta = atan2(Z-Z_center, Y-Y_center);

% keep unique values 
[Theta,ia,ib] = unique(Theta);
R = R(ia);
T = T(ia);
Rmean = mean(R);

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

% let's do some processing
Rs = smoothdata(R,'gaussian',2500);
[upperEnv,lowerEnv] = envelope(R, 100, 'peaks');

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rs,'k','linewidth',3);
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
plot(Theta,upperEnv,'r','linewidth',3);
plot(Theta,lowerEnv,'g','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

% remove outliers
a = 2; %  the parameter that define the threshold line :   th = (a*upperEnv+lowerEnv)/(1+a); % weigthed average curve using top and bottom envelopes

iter = 2; % nb of iterations 

for k = 1:iter
    th = (a*upperEnv+lowerEnv)/(1+a); % weigthed average curve using top and bottom envelopes
    % remove data below threshold th
    ind = R<=th;
    Theta(ind) = [];
    R(ind) = [];
    T(ind) = [];
    [upperEnv,lowerEnv] = envelope(R, 100, 'peaks'); % re -generate the envelopes for next iteration

end

% do a sine fit 
[Thetafit,Rfit,sol] = onesinefit(Theta,R);
% display sinus fit parameters in command window 
amplitude = sol(1)
frequency_Hz = sol(2)
phase_offset_rad = sol(3)
DC_offset = sol(4)

figure, scatter(Theta,R,15,T,'filled');colorbar
xlabel('Theta');
ylabel('R');
hold on
plot(Thetafit,Rfit,'--k','linewidth',2);
hold off

% final plot with interpolated T
Tfit = interp1(Theta,T,Thetafit);
figure, scatter(Thetafit,Rfit,15,Tfit,'filled');colorbar
xlabel('Theta');
ylabel('R');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tt,yp,sol] = onesinefit(t,y)
    % fit a sinewave 
    % the input data must contain only one period so there is only one min and one max
    % adapt otherwise
    %%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%
    ym = mean(y);                                             % Estimate offset
    [yu,iu] = max(y);
    [yl,il] = min(y);
    yr = (yu-yl)/2;                                             % Estimate Amplitude of ‘y’
    yz = y-ym;

    per = abs(t(iu) - t(il))*2;                                   % Estimate period
    fre = 1/per;                                           % Estimate FREQUENCY
    phe = -pi*(t(iu) + t(il))/per ;                                  % Estimate phase

    % stationnary sinus fit
    fit = @(b,x)  b(1) .* (sin(2*pi*x*b(2) + b(3))) + b(4);   % Objective Function to fit
    fcn = @(b) norm(fit(b,t) - y);                                                  % Least-Squares cost function
    sol = fminsearch(fcn, [yr;  fre;  phe;  ym;]);                                   % Minimise Least-Squares
    sol(3) = mod(sol(3),2*pi);
    tt = linspace(min(t),max(t),100);
    yp = fit(sol,tt);

end

[Ahmed]

Wow I didn't think of the sin wave method. I do see here that you removed tghe x coordiantes as well as took it to be a one layer thick cylinder (because it is with the data I have :p). My aim right now is to have a program that takes a multi layer thick cylinder and flatten out the inner and outer surfaces on a 2d plane (like you did). Your first picture made me realize my biggest error. The data is only on one surface so of course the inner and outer surfaces won't display properly since I'm ripping one surface into 2. Thank you very much for the code and the help, I'll request data that is actually multiple layers thick

[Mathieu NOE]

a better code with an improved sine fit function, that now works for sine wave having less or more than one period

requires a findpeak alike function : opted for peakseek which is much faster !!

attached the peakseek function or download it from the Fex :

PeakSeek - File Exchange - MATLAB Central

coords = readtable("Al_Case1.xlsx");
coords = table2array(coords);
X = coords(:,1);
Y = coords(:,2);
Z = coords(:,3);
T = coords(:,4);
Z_center = mean(Z);
Y_center = mean(Y);
R = ((Z-Z_center).^2 + (Y-Y_center).^2).^0.5;
Theta = atan2(Z-Z_center, Y-Y_center);

% keep unique values 
[Theta,ia,ib] = unique(Theta);
R = R(ia);
T = T(ia);
Rmean = mean(R);

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

% let's do some processing
Rs = smoothdata(R,'gaussian',2500);
[upperEnv,lowerEnv] = envelope(R, 100, 'peaks');

figure, scatter(Theta,R,15,T,'filled');colorbar
hold on
plot(Theta,Rs,'k','linewidth',3);
plot(Theta,Rmean*ones(size(Theta)),'k--','linewidth',3);
plot(Theta,upperEnv,'r','linewidth',3);
plot(Theta,lowerEnv,'g','linewidth',3);
xlabel('Theta');
ylabel('R');
hold off

% remove outliers
a = 2; %  the parameter that define the threshold line :   th = (a*upperEnv+lowerEnv)/(1+a); % weigthed average curve using top and bottom envelopes

iter = 2; % nb of iterations 

for k = 1:iter
    th = (a*upperEnv+lowerEnv)/(1+a); % weigthed average curve using top and bottom envelopes
    % remove data below threshold th
    ind = R<=th;
    Theta(ind) = [];
    R(ind) = [];
    T(ind) = [];
    [upperEnv,lowerEnv] = envelope(R, 100, 'peaks'); % re -generate the envelopes for next iteration

end

% do a sine fit 
[Thetafit,Rfit,sol] = sinefit(Theta,R);
% display sinus fit parameters in command window 
amplitude = sol(1)
frequency_Hz = sol(2)
phase_offset_rad = sol(3)
DC_offset = sol(4)

figure, scatter(Theta,R,15,T,'filled');colorbar
xlabel('Theta');
ylabel('R');
hold on
plot(Thetafit,Rfit,'--k','linewidth',2);
hold off

% final plot with interpolated T
Tfit = interp1(Theta,T,Thetafit);
figure, scatter(Thetafit,Rfit,15,Tfit,'filled');colorbar
xlabel('Theta');
ylabel('R');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tt,yp,sol] = sinefit(t,y)
    % fit a sinewave 
    % the input data can have a number of periods above or below one
    %%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%
    ym = mean(y);                                             % Estimate offset
    yz = y-ym;
    %yz = smoothdata(yz,'gaussian',10); % optionnal / adpat / remove
    [yu,iu] = max(y);
    [locsu, ~]=peakseek(yz,3,max(yz)/2);
    if isempty(locsu)
        locsu = iu;
    end
    
    [yl,il] = min(y);
    [locsl, ~]=peakseek(-yz,3,max(-yz)/2);
     if isempty(locsl)
        locsl = il;
     end
    
    til = t(locsl(1));
    tiu = t(locsu(1));
    
    yr = (yu-yl)/2;                                             % Estimate Amplitude 
    per = abs(til - tiu)*2;                                   % Estimate period
    fre = 1/per;                                           % Estimate FREQUENCY
    phe = pi-pi*(tiu + til)/per ;                                  % Estimate phase

    % stationnary sinus fit
    fit = @(b,x)  b(1) .* (sin(2*pi*x*b(2) + b(3))) + b(4);   % Objective Function to fit
    fcn = @(b) norm(fit(b,t) - y);                                                  % Least-Squares cost function
    sol = fminsearch(fcn, [yr;  fre;  phe;  ym;]);                                   % Minimise Least-Squares
    sol(3) = mod(sol(3),2*pi);
    tt = linspace(min(t),max(t),100);
    yp = fit(sol,tt);

end