I would like to draw a 4d isosurfce plot, where the 4th dimension, is the density of each grid voxel, and I would like it to be thee isovalue

M is a 3d grid of points.

Ro(t) is the density at grid point t.

Ro3D is thee density at each pixel as a 3D matrix that corrwesponds to each voxel.

How would you plot it?

Thank you!

M=DataSRX(find(idx==1),4:6);
StepSize=100; % [nm]
Xgrid=min(M(:,1)):StepSize:max(M(:,1));
Ygrid=min(M(:,2)):StepSize:max(M(:,2));
Zgrid=min(M(:,3)):StepSize:max(M(:,3));
Sigma=StepSize;
XGsize=numel(Xgrid); YGsize=numel(Ygrid); ZGsize=numel(Zgrid);
Z=zeros(1,XGsize*YGsize*ZGsize);
Ind=zeros(size(M,1),1);
Ro=Z;
% Z(1)=Xgrid(1),Ygrid(1),Zgrid(1); Z(2)=Xgrid(2),Ygrid(2),Zgrid(2);
% Z(Xn)=Xgrid(Xn),Ygrid(1),Zgrid(1); Z(Xn+1)=Xgrid(1),Ygrid(2),Zgrid(1);
Norm=1./((Sigma*sqrt(2*pi)^3));
for i=1:size(M,1)
    tx=ceil((M(i,1)-min(M(:,1)))/StepSize); ty=ceil((M(i,2)-min(M(:,2)))/StepSize); tz=ceil((M(i,3)-min(M(:,3)))/StepSize);
    if ty==0
        ty=1;
    end
    if tx==0
        tx=1;
    end
    if tz==0
        tz=1;
    end % if tz==0
    t=(XGsize+YGsize)*(tz-1)+XGsize*(ty-1)+tx;
    Z(t)=Z(t)+1;
    Ind(i)=t;
    Ro(t)=Norm*exp(((M(i,1)-Xgrid(tx))^2+(M(i,2)-Ygrid(ty))^2+(M(i,3)-Zgrid(tz))^2)/(2*Sigma^2))+Ro(t);
end % for i=1:size(M,1)

RoDense=Ro.*Z;
Ro3D=reshape(RoDense,[XGsize,YGsize,ZGsize]);
%{
Ro3D=zeros(XGsize,YGsize,ZGsize);
for i=1:numel(Ro)
   [I,J,K]=ind2sub([XGsize,YGsize,ZGsize],i);
   Ro3D(I,J,K)=Ro(i);
end % for i=1:numel(Ro)
%}
Ro3D_round=round(Ro3D,3);

[KSSV]

Not clear.......attaching data or an pictorial example with your expectations will help us to help you.

[Guy Nir]

Is the added code helpful to solving this question?

Thanks,

Guy

Answer 1

[Guy Nir]

To make it more clear I attached the data, and added code to parse the file.

% Plot density
%% Parse file
FileName='/Users/guynir/Dropbox (HMS)/Collaborations/MAMR Lab/Irene/CSV from paper/CSV_6columns/2017-08-03-19-34-25_Location-05_6column.csv';
[Header, DataSRX] = ReadSRXParticlesCSV(FileName);
figure(); scatter3(DataSRX(:,4),DataSRX(:,5),DataSRX(:,6),20,'filled')
axis equal
%% Separate homologs
Ploidy=2;
idx = kmeans(DataSRX(:,4:6),Ploidy);
figure(); scatter3(DataSRX(find(idx==1),4),DataSRX(find(idx==1),5),DataSRX(find(idx==1),6),20,'filled','r')
hold on
scatter3(DataSRX(find(idx==2),4),DataSRX(find(idx==2),5),DataSRX(find(idx==2),6),20,'filled','b')
axis equal
%% Surface plots
M=DataSRX(find(idx==1),4:6);
shp = alphaShape(M);
figure(); plot(shp)

%% Density map

StepSize=100; % [nm]
Xgrid=min(M(:,1)):StepSize:max(M(:,1));
Ygrid=min(M(:,2)):StepSize:max(M(:,2));
Zgrid=min(M(:,3)):StepSize:max(M(:,3));
Sigma=StepSize;
XGsize=numel(Xgrid); YGsize=numel(Ygrid); ZGsize=numel(Zgrid);
Z=zeros(1,XGsize*YGsize*ZGsize);
Ind=zeros(size(M,1),1);
Ro=Z;
% Z(1)=Xgrid(1),Ygrid(1),Zgrid(1); Z(2)=Xgrid(2),Ygrid(2),Zgrid(2);
% Z(Xn)=Xgrid(Xn),Ygrid(1),Zgrid(1); Z(Xn+1)=Xgrid(1),Ygrid(2),Zgrid(1);
Norm=1./((Sigma*sqrt(2*pi)^3));
for i=1:size(M,1)
    tx=ceil((M(i,1)-min(M(:,1)))/StepSize); ty=ceil((M(i,2)-min(M(:,2)))/StepSize); tz=ceil((M(i,3)-min(M(:,3)))/StepSize);
    if ty==0
        ty=1;
    end
    if tx==0
        tx=1;
    end
    if tz==0
        tz=1;
    end % if tz==0
    t=(XGsize+YGsize)*(tz-1)+XGsize*(ty-1)+tx;
    Z(t)=Z(t)+1;
    Ind(i)=t;
    Ro(t)=Norm*exp(((M(i,1)-Xgrid(tx))^2+(M(i,2)-Ygrid(ty))^2+(M(i,3)-Zgrid(tz))^2)/(2*Sigma^2))+Ro(t);
end % for i=1:size(M,1)

RoDense=Ro.*Z;
Ro3D=reshape(RoDense,[XGsize,YGsize,ZGsize]);
Ro3D_round=round(Ro3D,3);