3D Histogram from raw data

16 views (last 30 days)
Biraj Khanal
Biraj Khanal on 17 May 2023
Commented: Biraj Khanal on 24 May 2023
I am trying to create the a rainflow matrix out of time history of a load data. I can do it easily with the rainflow command but I am trying to learn how I can produce it myself. It is basically a 3D Histogram of three variables ( Mean, Range and Cycle count , random X,Y and Z let's say )
This is what a rainflow matrix looks like for reference.
For any dataset with X,Y and Z data, I undertand that the first step is to discretize X and Y . What I do not understand is how I can allot the Z value to each bins and ultimately plot something like this.

Answers (1)

Mathieu NOE
Mathieu NOE on 17 May 2023
hello
see example below - this is equivalent to hist3
% dummy data (col vectors)
x = rand(100,1);
y = rand(100,1);
z = rand(100,1);
%% bin the data (Hist3 clone)
nBins = 10; %number of bins (there will be nBins + 1 edges)
XEdge = linspace(min(x),max(x),nBins+1);
YEdge = linspace(min(y),max(y),nBins+1);
[~, xBin] = histc(x, XEdge);
[~, yBin] = histc(y, YEdge);
% count number of elements per (x,y) pair
[xIdx, yIdx] = meshgrid(1:nBins, 1:nBins);
xyPairs = [xIdx(:), yIdx(:)];
Z = zeros(size(xyPairs,1),1);
for i = 1:size(xyPairs,1)
Z(i) = sum(ismember([xBin, yBin], xyPairs(i,:), 'rows'));
end
% Reshape nPerBin to grid size
Z = reshape(Z, [nBins, nBins]);
% plot data
figure(1)
b = bar3(Z);
% change bar color to match Z value
for k = 1:length(b)
zdata = b(k).ZData;
b(k).CData = zdata;
b(k).FaceColor = 'interp';
end
colorbar('vert');
% Change colormap
colormap jet
% Label the axes
xlabel('x')
ylabel('y')
title('hist3 simulation');
  2 Comments
Biraj Khanal
Biraj Khanal on 19 May 2023
Edited: Biraj Khanal on 19 May 2023
Hi Mathieu , this was really helpful to get started. I replaced the histc with histcounts ( which I believe is the same thing) . However, what you have done is basically counted the xy pairs in each bins and presented Z axis as the count. Note that you have not used the initial z data at all. What I am doing is trying to put those respective z data (because the data is a x,y,z pair) in each bin and make a sum of them to present as Z of the final histogram. I have done that changing a line of your code. Check the for loop.
%dummy data (col vectors)
x = rand(100,1);
y = rand(100,1);
z = rand(100,1);
%% bin the data (Hist3 clone)
nBins = 10; %number of bins (there will be nBins + 1 edges)
% XEdge = linspace(min(x),max(x),nBins+1);
% YEdge = linspace(min(y),max(y),nBins+1);
% [~, xBin] = histc(x, XEdge);
% [~, yBin] = histc(y, YEdge);
[~,XEdge,XBin] = histcounts(x,nBins);
[~,YEdge,YBin] = histcounts(y,nBins);
[xIdx, yIdx] = meshgrid(1:nBins, 1:nBins);
xyPairs = [xIdx(:), yIdx(:)];
Z = zeros(size(xyPairs,1),1);
for i = 1:size(xyPairs,1)
%Z(i) = sum(ismember([xBin, yBin], xyPairs(i,:), 'rows'));
z_idx = find(ismember([XBin, YBin], xyPairs(i,:), 'rows'));
Z(i) = sum(z(z_idx));
end
Z = reshape(Z, [nBins, nBins]);
%Z = flip(Z,1);
% plot data
figure(1)
b = bar3(Z);
% change bar color to match Z value
for k = 1:length(b)
zdata = b(k).ZData;
b(k).CData = zdata;
b(k).FaceColor = 'interp';
end
colorbar('vert');
% Change colormap
colormap jet
% Label the axes
xlabel('x')
ylabel('y')
title('hist3 simulation');
Now, there are couple of more things that I am trying to do:
  1. As you can see, the mesh in XY is not as we generally liike to see, The Y bins should be flipped. If only for the sake of Z, I could just flip Z. But I believe there is a better way to arrange this.
  2. Instead of the bin index, I want to see the edges so that the histogram is relatable to the original data.
Do you have some suggestion?
Biraj Khanal
Biraj Khanal on 24 May 2023
For anyone having similar problem, here is what I did to address the issues in my comment above :
Z = flip(Z,1); %flips the bar graph
ylim([0 nBins+1 ]); % to make the grid well viyible
xlim([0 nBins+1]);
tick_num = 5; %set the number of ticks on x and y axis
x = 0: nBins/tick_num : nBins ; x = x+0.4 ; x(1) =0.6; %first bin is plotted from x = 0.6 , each bin ends at edge+0.4
y = x;
xticks(x)
yticks(y)
xticklabels([XEdge(1: nBins/tick_num :end)]);
yticklabels(flip([YEdge(1: nBins/tick_num :end)]));
If anyone can suggest simpler way to go about this, it would be great!!

Sign in to comment.

Categories

Find more on Data Distribution Plots in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!