Added distribution sampling when direct pixel replication would produce over a million points.

srcx/MATLAB/PixelReplicate.m

@@ -75,3 +75,4 @@ ind = unique(ind);
 bw(ind) = true;
 bwd = bwdist(~bw);
 end
+

srcx/MATLAB/PixelReplicateSample.m

+% ptsReplicate_xy = PixelReplicateSample(im,pct, label)
+%
+% im can be a logical image - if so it should contain only one connected
+%   component
+% im can be a label image - it can contain any number of connected components 
+%   and label is used to select the component to replicate
+%
+% pct is the percentage reduction of points to use relative to the total
+%   distance pixel replication generates
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Copyright (c) 2016, Drexel University
+% All rights reserved.
+% 
+% Redistribution and use in source and binary forms, with or without
+% modification, are permitted provided that the following conditions are met:
+% 
+% * Redistributions of source code must retain the above copyright notice, this
+%   list of conditions and the following disclaimer.
+% 
+% * Redistributions in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the documentation
+%   and/or other materials provided with the distribution.
+% 
+% * Neither the name of PixelRep nor the names of its
+%   contributors may be used to endorse or promote products derived from
+%   this software without specific prior written permission.
+% 
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function ptsReplicate_xy = PixelReplicateSample(im,pct, label)
+    if ( ~exist('pct','var') )
+        pct = 0.2;
+    end
+
+    if ( ~exist('label','var') )
+        label = [];
+    end
+
+    bw = (im > 0);
+    if ( ~isempty(label) )
+        bw = (im==label);
+    end 
+
+    bwd = bwdist(~bw);
+    npts = sum(bwd(:));
+    
+    ptsReplicate_xy = RejectionSim(round(pct*npts), bwd);
+end

srcx/MATLAB/RejectionSim.m

+% X = RejectionSim(n,pdfIm)
+% 
+% Simulate n points from the discrete empirical density pdfIm using the
+% rejection method.
+
+function X = RejectionSim(n,pdfIm)
+    d = ndims(pdfIm);
+    
+    % Generate only near non-zero pixels
+    bwIm = imdilate(pdfIm > 0, strel('square',3));
+    validIdx = find(bwIm);
+    
+    coordCell = cell(1,d);
+    [coordCell{:}] = ind2sub(size(pdfIm), validIdx);
+    validCoords = [coordCell{:}];
+    
+    % Guarantee pdfIm is distribution and find probability of rejection c
+    
+    % g is the pdf of a uniform distribution over the valid pixels
+    g = 1/length(validIdx);
+    pdfIm = pdfIm / sum(pdfIm(:));
+    
+    c = max(pdfIm(:)) / g;
+    
+    lastIdx = 0;
+    X = zeros(n,d);
+    for i=1:ceil(2*c)
+        U = rand(n,1);
+        coordIdx = randi(length(validIdx),n,1);
+        C = validCoords(coordIdx,:) + rand(n,d) - 0.5;
+        
+        cellCoord = num2cell(C,1);
+        f = interpn(pdfIm, cellCoord{:}, 'linear');
+        
+        bKeep = (U <= f/(c*g));
+        
+        numValid = nnz(bKeep);
+        writeNum = min(lastIdx+numValid,n)-lastIdx;
+        
+        validC = C(bKeep,:);
+        X(lastIdx+(1:writeNum),:) = validC(1:writeNum,:);
+        
+        lastIdx = lastIdx + writeNum;
+        if ( lastIdx >= n )
+            break;
+        end
+    end
+    
+    if ( lastIdx < n )
+        warning(['Only generated ' num2str(lastIdx) ' valid points!']);
+    end
+    
+    swapCoord = [2 1 3:size(X,2)];
+    X = X(:,swapCoord);
+end

srcx/MATLAB/SurfHist.m

+function binIm = SurfHist(X, binSz)
+    if ( ~exist('binSz','var') )
+        binSz = 1;
+    end
+    
+    if ( isscalar(binSz) )
+        binSz = repmat(binSz,1,2);
+    end
+
+    dataRange = [min(X,[],1); max(X,[],1)];
+    dataSize = diff(dataRange,1) + binSz;
+    
+    totalSize = ceil(dataSize ./ binSz);
+    binIm = zeros(totalSize([2 1]));
+    
+    shiftData = X-repmat(dataRange(1,:)-binSz,size(X,1),1);
+    binData = round(shiftData./repmat(binSz,size(X,1),1));
+    
+    [uqBins,~,ic] = unique(binData(:,[2 1]),'rows');
+    binCell = num2cell(uqBins,1);
+    
+    uqIdx = sub2ind(size(binIm), binCell{:});
+    
+    dataMap = sort(ic);
+    repPts = [1;diff(dataMap)];
+    uqCounts = diff([find(repPts);length(dataMap)+1]);
+    
+    binIm(uqIdx) = uqCounts;
+    
+    binRange = cell(1,length(binSz));
+    for i=1:length(binSz)
+        startBin = dataRange(1,i)-binSz(i)/2;
+        endBin = dataRange(2,i)+binSz(i)/2;
+        binRange{i} = startBin:binSz(i):(endBin-binSz(i)/2);
+    end
+    
+    [XX,YY] = meshgrid(binRange{:});
+    surf(XX,YY,binIm, 'EdgeColor','none');
+end

srcx/MATLAB/demoPR.m

@@ -41,6 +41,7 @@ cmap=hsv(K);
 % generate a random image. you can use your own image here instead...
 bw=GetRandomEllipseImage(K); 
 
+figure(1);
 hold off
 imshow(bw)
 RMAX=K*30;
@@ -57,7 +58,13 @@ hold on
 %   ....
 % end
 
+% For more than a million replicate points use sampling at 10% instead.
+numPoints = sum(round(bwdist(~bw)));
+if ( numPoints < 1e6  )
     ptsReplicated = PixelReplicate(bw);
+else
+    ptsReplicated = PixelReplicateSample(bw,0.1);
+end
 
 % fit gmm to PR points
 % NOTE - more replicates is more accurate fit, but takes longer. you can
@@ -66,10 +73,10 @@ objPR = fitgmdist(ptsReplicated, K, 'replicates',5);
 
 % draw the clustering
 [y,x]=find(bw);
-idx = objPR.cluster([y,x]);
+idx = objPR.cluster([x,y]);
 for kk=1:K
     idxK=find(idx==kk);
-    plot(x(idxK),y(idxK),'.')
+    plot(x(idxK),y(idxK),'.','color',cmap(kk,:))
 end
 
 %draw the gaussian isocontours for the triangular elliptical approximation