如何使用给定数据绘制 3D 数密度散点图?

how to do 3D number-density scatter plot using given data?

我有一个立方体盒子,每个方向的大小为 300。我有一些数据(包含 X、Y、Z 坐标)代表此框中近 100 万个数据点的分布。我想为它们的数字密度指定一种颜色(在这种情况下,它是一个密集的数量,用于描述可数对象的集中程度数据点)。换句话说,使用颜色来说明哪个部分在数据点方面比其他部分更浓缩。最终图像中颜色条的索引应表示使用该颜色指定的数据点的百分比。

我试图通过将立方体中的整个 space 分成 100 万个更小的立方体(每个立方体在所有方向上的长度均为 3)来做到这一点。通过计算这些立方体中的粒子数量,我将知道它们在盒子中的分布情况以及其中存在的数据点数量。然后我可以为他们指定一种我没有成功计算和指定的颜色。任何建议表示赞赏。

%reading the files
[FileName,PathName,FilterIndex] = uigetfile('H:\*.txt','MultiSelect','on');
numfiles = size(FileName,2);%('C:\final1.txt');
j=1;
X=linspace(0,300,100);
for ii = 1:numfiles
    FileName{ii}
    entirefile = fullfile(PathName,FileName{ii});
    a = importdata(entirefile);
    x = a(:,2);
    y = a(:,3);
    z = a(:,4);
    %% I DON'T KNOW HOW TO CREAT THIS LOOP TO COUNT FOR THE NUMBER OF PARTICLES WITHIN EACH DEFINED CUBE %%
    for jj = 2:size(X,2) 
        %for kk=1:m
        if x(:)<X(jj) & y(:)<X(jj) & z(:)<X(jj)
           x;
        end
        %end
    end
    h=figure(j);
    scatter3(x, y, z, 'filled', 'MarkerSize', 20);
    cb = colorbar();
    cb.Label.String = 'Probability density estimate';
end

我需要得到类似下图的结果。但我需要每种颜色指定的数据点百分比。提前致谢。

这里是link一个采样数据。

这是一种计算点云 3D 密度的方法。尽管我担心您提供的示例数据不会产生与您的示例图像相同的 3D 分布。

为了计算密度,方法分为几个步骤:

  • 计算 [X,Y] 平面中的二维密度:计算每个 (x,y) bin 中的点数。但是,在这个阶段,这个点数包含给定 bin 的所有 Z 列。
  • 对于每个 非空 (x,y) bin,计算沿 Z 列的分布。我们现在有落在每个 (x,y,z) bin 中的点数。计算 density/percentage 只需将每个计数除以总点数即可。
  • 现在,对于每个 非空 (x,y,z) bin,我们确定属于该 bin 的点的线性索引。然后,我们将 bin 值(颜色、百分比或与此 bin 关联的任何值)分配给所有已识别的点。
  • 显示结果。

在代码中,它是这样的:

%% Import sample data
entirefile = '1565015520323.txt' ;
a = importdata(entirefile);
x = a(:,1);
y = a(:,2);
z = a(:,3);
npt = numel(x) ; % Total Number of Points

%% Define domain and grid parameters
nbins    = 100 ;
maxDim   = 300 ;
binEdges = linspace(0,maxDim,nbins+1) ;

%% Count density
% we start counting density along in the [X,Y] plane (Z axis aglomerated)
[Nz,binEdges,~,binX,binY] = histcounts2(y,x,binEdges,binEdges) ;

% preallocate 3D containers
N3d = zeros(nbins,nbins,nbins) ; % 3D matrix containing the counts
Npc = zeros(nbins,nbins,nbins) ; % 3D matrix containing the percentages
colorpc = zeros(npt,1) ;         % 1D vector containing the percentages

% we do not want to loop on every block of the domain because:
%   - depending on the grid size there can be many
%   - a large number of them can be empty
% So we first find the [X,Y] blocks which are not empty, we'll only loop on
% these blocks.
validbins = find(Nz) ;                              % find the indices of non-empty blocks
[xbins,ybins] = ind2sub([nbins,nbins],validbins) ;  % convert linear indices to 2d indices
nv = numel(xbins) ;                                 % number of block to process

% Now for each [X,Y] block, we get the distribution over a [Z] column and
% assign the results to the full 3D matrices
for k=1:nv
    % this block coordinates
    xbin = xbins(k) ;
    ybin = ybins(k) ;

    % find linear indices of the `x` and `y` values which are located into this block
    idx = find( binX==xbin & binY==ybin ) ;
    % make a subset with the corresponding 'z' value
    subZ = z(idx) ;
    % find the distribution and assign to 3D matrices
    [Nz,~,zbins] = histcounts( subZ , binEdges ) ;
    N3d(xbin,ybin,:) = Nz ;         % total counts for this block
    Npc(xbin,ybin,:) = Nz ./ npt ;  % density % for this block

    % Now we have to assign this value (color or percentage) to all the points
    % which were found in the blocks
    vzbins = find(Nz) ;
    for kz=1:numel(vzbins)
        thisColorpc = Nz(vzbins(kz)) ./ npt * 100 ;
        idz   = find( zbins==vzbins(kz) ) ;
        idx3d = idx(idz) ;
        colorpc(idx3d) = thisColorpc ;
    end

end
assert( sum(sum(sum(N3d))) == npt ) % double check we counted everything

%% Display final result
h=figure;
hs=scatter3(x, y, z, 3 , colorpc ,'filled' );
xlabel('X'),ylabel('Y'),zlabel('Z')
cb = colorbar ;
cb.Label.String = 'Probability density estimate';

正如我一开始所说,结果与您的示例图像略有不同。此样本集产生以下分布:

如果你想 "double check" 结果不是垃圾,你可以查看每个轴上的二维密度结果,并检查它是否与你的点的表观分布相匹配:

%% Verify on 3 axis:
Nz = histcounts2(y,x,binEdges,binEdges) ./ npt *100 ;
Nx = histcounts2(z,y,binEdges,binEdges) ./ npt *100 ;
Ny = histcounts2(x,z,binEdges,binEdges) ./ npt *100 ;
figure
ax1=subplot(1,3,1) ; bz = plotDensity(Nz,ax1) ; xlabel('X'),ylabel('Y') ;
ax2=subplot(1,3,2) ; bx = plotDensity(Nx,ax2) ; xlabel('Y'),ylabel('Z') ;
ax3=subplot(1,3,3) ; by = plotDensity(Ny,ax3) ; xlabel('Z'),ylabel('X') ;

点击图片放大:

plotDensity.m的代码:

function hp = plotDensity(Ndist,hax)
    if nargin<2 ; hax = axes ; end
    hp = bar3(Ndist,'Parent',hax) ;
    for k = 1:length(hp)
        zdata = hp(k).ZData;
        hp(k).CData = zdata;
        hp(k).FaceColor = 'interp';
    end
    shading interp