cleaned up phantom code

src/MATLAB/phantomKymograph/createPhantom.m

+function phantomKymo = createPhantom(phSize,phTracks,phVel_mu,phVel_sigma)
+
+tmp = zeros(phSize,'uint8');
+
+x = randsample(1:25,phTracks);
+y = randsample(1:25,phTracks);
+
+z = ones(phTracks,1);
+trackSrcPts = [x', y', z];
+
+for i = 1:length(trackSrcPts)
+    
+    b = trackSrcPts(i,:);
+    
+    phVelocity = abs(normrnd(phVel_mu,phVel_sigma));
+
+    % if round(phVelocity) == 0
+        % phVelocity = 1;
+    % end
+
+    rVal = randi(255);
+    tmp(b(1),b(2),b(3)) = phVelocity;
+    tmp(b(1),b(2),b(3),2) = rVal;
+
+    for j = trackSrcPts(i,3) + 1:phSize(3)
+
+        X = b + round(phVelocity);
+
+        if any(X > phSize(1))
+            break
+        end
+
+        tmp(X(1), X(2), j) = phVelocity;
+        tmp(X(1), X(2), j, 2) = rVal;
+
+        b = X;
+
+    end
+
+end
+
+phantomKymo = tmp;
+
+end
\ No newline at end of file

src/MATLAB/phantomKymograph/createSyntheticDataset.m

+function syntheticDataset = createSyntheticDataset(datasetParameters)
+
+syntheticDataset = [];
+
+phSize = datasetParameters.phantomSize;
+phTracks = datasetParameters.phantomTracks;
+
+for i = 1:length(datasetParameters.phantomVelocity_mu)
+
+    phVel_mu = datasetParameters.phantomVelocity_mu(i);
+    phVel_sigma = datasetParameters.phantomVelocity_sigma(i);
+
+    for j = 1:datasetParameters.numKymos
+
+        tmp = [];
+        tmp.phantomName = strcat(num2str(i),'_',num2str(j));
+        tmp.phantomGroup = i;
+        tmp.datetime = string(datetime('now'));
+        tmp.phantomSize = phSize;
+        tmp.numTracks = phTracks;
+        tmp.phantomVelocity = [phVel_mu;phVel_sigma];
+        tmp.phantomKymo = createPhantom(phSize,phTracks,phVel_mu,phVel_sigma);
+
+        syntheticDataset = [syntheticDataset;tmp];
+
+    end
+end
+
+end
\ No newline at end of file

src/MATLAB/phantomKymograph/drawSyntheticSSF.m

+function drawSyntheticSSF(A,K,syntheticDataset,plotTitle)
+
+[~,Y] = Cluster.SpectralCluster(A,K);
+
+[~,~,idxTrue] = unique([syntheticDataset.phantomGroup]');
+[mm,ss] = CSF.csf_spatial(Y,idxTrue);
+
+sym = {'o','x','*'};
+
+hold on
+hx = [];
+
+for ff=1:size(Y,1)
+    if syntheticDataset(ff).phantomGroup == 1
+        color = [1 0 0];
+        idxSym = 1;
+    elseif syntheticDataset(ff).phantomGroup == 2
+        color = [0 1 0];
+        idxSym = 2;
+    else
+        color = [0 0.4470 0.7410];
+        idxSym = 3;
+
+    end
+
+    % h1 = plot3(Y(ff,1),Y(ff,2),Y(ff,3),'.','color',color,'UserData',ff,'MarkerSize',8); % for debugging, does not use different markers
+    h1 = plot3(Y(ff,1),Y(ff,2),Y(ff,3),sym{idxSym},'color',color,'MarkerSize',6);
+    hx(idxSym) = h1;
+
+    % for debugging, lists names for each point
+    % phantomID = strcat({' '},syntheticDataset(ff).phantomName);
+    % text(Y(ff,1),Y(ff,2),Y(ff,3),phantomID,'Interpreter','none','color',color, 'FontSize',8)
+
+end
+
+axis square
+
+set(gcf,'color','w')
+xlabel('k1');ylabel('k2'),zlabel('k3')
+lgd = legend(hx,{'v=1','v=3','v=5'});
+set(lgd,'color','none');
+
+title([plotTitle,' = [', num2str(round(mm,2,'significant')), ',' ...
+    ,num2str(round(ss,2,'significant')),']']);
\ No newline at end of file

src/MATLAB/phantomKymograph/getPhantomClipLimits.m

+function clipLimits = getPhantomClipLimits(syntheticDataset,targetChannel)
+
+for ff = 1:length(syntheticDataset)
+    im_ff = syntheticDataset(ff).phantomKymo(:,:,:,targetChannel);
+    
+    kp = im_ff(find(im_ff));
+
+    if ff == 1
+        clipLimits(targetChannel,:) = [min(kp),max(kp)];
+    else
+        clipLimits(targetChannel,1) = min(clipLimits(targetChannel,1),min(kp));
+        clipLimits(targetChannel,2) = max(clipLimits(targetChannel,2),max(kp));
+    end
+end
+
+4;
\ No newline at end of file

src/MATLAB/phantomKymograph/goFigure.m

+% load('syntheticDataset_16.mat');
+
+K = 3;
+
+% tmp_xy(:,:,1) = squeeze(max(syntheticDataset(8).phantomKymo(:,:,:,1),[],2));
+% tmp_xy(:,:,2) = squeeze(max(syntheticDataset(51).phantomKymo(:,:,:,1),[],2));
+% tmp_xy(:,:,3) = squeeze(max(syntheticDataset(101).phantomKymo(:,:,:,1),[],2));
+
+tmp_xy(:,:,1) = squeeze(max(syntheticDataset(1).phantomKymo(:,:,:,1),[],2));
+tmp_xy(:,:,2) = squeeze(max(syntheticDataset(11).phantomKymo(:,:,:,1),[],2));
+tmp_xy(:,:,3) = squeeze(max(syntheticDataset(21).phantomKymo(:,:,:,1),[],2));
+
+figure('Position',[300,200,900,300]);
+layout1 = tiledlayout(1,3,'TileSpacing','compact','Padding','tight');
+layout2 = tiledlayout(layout1,1,3,'TileSpacing','tight','Padding','tight');
+layout2.Layout.Tile = 1;
+layout2.Layout.TileSpan = [1 1];
+
+axVel = nexttile(layout1,2,[1 1]);
+plotTitle = 'CSF_{vel}^{3}';
+drawSyntheticSSF(syntheticDatasetVel_d,K,syntheticDataset,plotTitle);
+view(axVel,-40,20)
+
+axRval = nexttile(layout1,3,[1 1]);
+plotTitle = 'CSF_{rnd}^{3}';
+drawSyntheticSSF(syntheticDatasetRVal_d,K,syntheticDataset,plotTitle);
+view(axRval,-40,20)
+
+for i = 1:3
+    ax = nexttile(layout2);
+    imagesc(tmp_xy(:,:,i));
+    ylim([1 160])
+    % ylim([1 180])
+    xlim([1 24])
+    yticks([30 80 130]);
+    xticks([6 12 18]);
+    clim([0,7]);
+
+    if i == 1
+        ax.YLabel.String = 'xy \Downarrow x';
+    end
+
+    if i == 2
+        ax.XLabel.String = 'time';
+    end
+
+end
+
+cm = colormap();    
+cm(1,:) = [1,1,1];
+colormap(cm);
+cb = colorbar;
+cb.Label.String = 'Velocity SSF';
+
+% outfile = 'phantom.png';
+% exportgraphics(gcf,outfile,'Resolution',300);
\ No newline at end of file

src/MATLAB/phantomKymograph/goSyntheticDataset.m

+clear;
+
+% velocity mean and ss determine how large the phantom must be
+
+synthDatasetParameters.phantomSize = [200,200,25];
+synthDatasetParameters.phantomTracks = 10;
+synthDatasetParameters.phantomVelocity_mu = [1,3,5];
+synthDatasetParameters.phantomVelocity_sigma = [.5,.5,.5];
+synthDatasetParameters.numKymos = 50;
+
+syntheticDataset = createSyntheticDataset(synthDatasetParameters);
+
+syntheticDatasetVel_d = goSyntheticSSF(syntheticDataset,1);
+syntheticDatasetRVal_d = goSyntheticSSF(syntheticDataset,2);
+
+save('syntheticDataset_01.mat','-v7.3');
\ No newline at end of file

src/MATLAB/phantomKymograph/goSyntheticSSF.m

+function d = goSyntheticSSF(syntheticDataset,targetChannel)
+
+tStart = tic();
+
+d = zeros(length(syntheticDataset));
+p = ljsStartParallel();
+
+cmdList = NCD.dParallelCommandList(ones(length(syntheticDataset)),p.NumWorkers);
+W = size(cmdList,2); % for slicing j in parfor
+H = size(cmdList,1);
+dxx = zeros(H,W);
+parfor i=1:H
+    
+    for j = 1:W
+
+        t0 = tic();
+        ff = cmdList(i,j,2); % ff is the row var -- repeats
+        gg = cmdList(i,j,1);
+
+        if 0 == ff || 0 == gg
+            continue
+        end
+
+        im_ff = syntheticDataset(ff).phantomKymo(:,:,:,targetChannel);
+
+        if ff == gg
+            d1 = SSF.ncd_ssf_volume(im_ff,im_ff);
+            dxx(i,j) = d1;
+        else
+            im_ff = syntheticDataset(ff).phantomKymo(:,:,:,targetChannel);
+            im_gg = syntheticDataset(gg).phantomKymo(:,:,:,targetChannel);
+
+            d1 = SSF.ncd_ssf_volume(im_ff,im_gg);
+            d2 = SSF.ncd_ssf_volume(im_gg,im_ff);
+
+            dxx(i,j) = min(d1,d2);
+        end
+        tElapsed = toc(t0);
+        fprintf('%s tElapsed = %0.2f\n',mat2str([ff,gg,dxx(i,j)]),tElapsed);
+
+    end
+end
+
+% convert dxx to distance matrix d
+for i = 1 : size(dxx,1)
+    for j = 1 : size(dxx,2)
+        c = cmdList(i,j,1);
+        r = cmdList(i,j,2);
+        if 0 == r || 0 == c
+            continue
+        end
+        d(r,c) = dxx(i,j);
+    end
+end
+
+for i = 1 : size(d,1)
+    for j = 1 : size(d,2)
+        if 0 == d(i,j)
+            d(i,j) = d(j,i);
+        end
+    end
+end
+tElapsed = toc(tStart);
+fprintf('total time tElapsed = %0.2f\n',tElapsed);
+
+end
\ No newline at end of file

src/MATLAB/phantomSSF/ssfPhantom.m

+function [cnRatio,ssfValue] = ssfPhantom(density)
+
+% set width/length of the image
+wl = 200;
+
+imRef = zeros(wl,wl);
+
+xc = wl/2;
+yc = wl/2;
+
+r = wl/10;
+
+for x = 1:size(imRef,1)
+    for y = 1:size(imRef,2)
+ 
+        d = sqrt((x-xc)^2+(y-yc)^2);
+        if d >= r
+            imRef(x,y) = 1;
+        end
+        
+    end
+end
+
+% imRef = imdilate(imRef,strel('disk',2));
+
+bw1 = bwdist(imRef);
+bw2 = bwdist(~imRef);
+
+% calculate C/N
+idxNuclear = find(bw1 > 0 & bw1 < 3);
+idxCyto = find(bw2 > 0 & bw2 < 3);
+
+imNoise = imnoise(imRef,'salt & pepper',density);
+
+% length(find(imRef ~= imNoise))/length(imNoise(:))
+
+% imNoise = imdilate(imRef,strel('disk',density*10));
+% imNoise = imerode(imRef,strel('disk',density*10));
+
+logRadius = ((1/sqrt(2)).*r);
+logRadius = repelem(logRadius,3);
+
+targetDevice = getCudaTarget();
+imNoiseLoG = HIP.LoG(imNoise(:,:,1),logRadius,targetDevice);
+
+% cnRatio = mean(imNoise(idxNuclear))/mean(imNoise(idxCyto));
+cnRatio = mean(imNoise(idxCyto))/mean(imNoise(idxNuclear));
+ssfValue = max(imNoiseLoG(:));
+
+end
\ No newline at end of file