Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % generate data
- tempdisk = strel('disk',922/2);
- tempdisk = double(repmat((1+sqrt(-1)).*tempdisk.Neighborhood,[1 1 15]));
- tempnoise = (rand(921,921,15)+sqrt(-1).*rand(921,921,15))./10;
- tempim1 = double(imresize(mean(imread('cameraman.tif'),3),[921,921]));
- tempim1 = repmat(tempim1./max(tempim1(:)),[1 1 15]);
- tempim2 = double(rgb2hsv(imread('fabric.png')));
- tempim2 = imresize(tempim2(:,:,1),[921,921]);
- tempim2 = repmat(tempim2./max(tempim2(:)),[1 1 15]);
- sin1 = repmat(permute(sin(2.*pi.*(0:14)./15),[1 3 2]),[921 921 1]);
- cdat = (sin1.*tempim1.*exp(-sqrt(-1).*2.*pi.*sin1.*(tempim2.^2)).*tempdisk+tempnoise);
- %this is what the mean data looks like
- meanm = mean(abs(cdat),3);
- meanp = angle(mean(cdat,3));
- figure; subplot(1,2,1); imshow(meanm,[]); title('mean magnitude'); subplot(1,2,2); imshow(meanp,[]); title('mean phase')
- %get all pixel vectors in a single matrix
- cvect = reshape(permute(cdat, [3,1,2]), [15, 921*921]);
- %k means and eigs dont accept complex, so convert to real here?
- cvectT = [real(cvect);imag(cvect)]';
- %lets say 10000 by 10000 matrices are still ok
- nvox = 10000;
- nslab = floor(length(cvectT)/nvox);
- nrest = rem(length(cvectT), nvox);
- % spectral clusetring according to http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf
- keig = 50;%how many eigenvectors needed? more is better
- affinity_sigma = 1;% i dont understand how to calculate this either
- tic
- for islab = (nslab+1):-1:1;
- islab/nslab
- toc
- if islab>nslab
- voxrange = (1:nrest) + ((islab-1)*nvox);;
- else
- voxrange = (1:nvox) + ((islab-1)*nvox);
- end
- Aff = exp( -squareform(pdist(cvectT(voxrange,:))).^2/(2*affinity_sigma^2) ); % affinity matrix
- Dsq = sparse(size(Aff,1),size(Aff,2)); %degree matrix
- for idiag=1:size(Aff,1)
- Dsq(idiag,idiag) = sum(Aff(idiag,:))^(1/2);
- end
- Lap = Dsq * Aff * Dsq; %normalize affinity matrix
- [eigVectors(voxrange,1:keig), eigValues] = eigs(Lap, keig); %eigenvectors of normalized aff mat
- normEigVectors(voxrange,1:keig) = eigVectors(voxrange,1:keig)./repmat(sqrt(sum(abs(eigVectors(voxrange,1:keig)).^2,2)), [1 keig]); %normalize rows of eigen vectors, normr only works on real numbers
- [idx,C,sumd,D] = kmeans([real(normEigVectors(voxrange,1:keig)),imag(normEigVectors(voxrange,1:keig))], 5); %k means on normalized eigenvecotrs
- idxval(voxrange) = idx;
- end
- idxim = reshape(idxval, [1280, 1280]);
- figure; imshow(idxim,[])
- toc
Add Comment
Please, Sign In to add comment