I'm currently trying to implement a method to generate TSP art, and for that I need a list of points (x,y), the local density of which is proportional to the gray scale pixel value of a given image.
My first thought was: well that works pretty much like Inverse Transform Sampling for statistics (you want to draw a sample that matches a given probability density function but you can only create a sample that is uniformly distributed).
I implemented this and it works fairly well, as evident by executing this code:
%% Load image, adjust it for our needs
im=imread('http://goo.gl/DDwV3t'); %load random headshot from google
im=imadjust(im,stretchlim(im,[.01,.65]),[]);
im=im2double(rgb2gray(im));
im=im(10:end-5,50:end-5);
figure;imshow(im);title('original');
im=1-im; %we want black dots on white background
im=flipud(im); %and we want it the right way up
%% process per row
imrow = cumsum(im,2);
imrow=imrow*size(imrow,1)./repmat(max(imrow,[],2),1,size(imrow,2));
y=1:size(imrow,2);
ximrow_i = zeros(size(imrow));
for i = 1:size(imrow,1)
mask =logical([diff(imrow(i,:))>=0.01,0]); %needed for interp
ximrow_i(i,:) = interp1(imrow(i,mask),y(mask),y);
end
y=1:size(ximrow_i,1);
y=repmat(y',1,size(ximrow_i,2));
y1=y(1:5:end,1:5:end); %downscale a bit
ximcol_i1=ximrow_i(1:5:end,1:5:end); %downscale a bit
figure('Color','w');plot(ximcol_i1(:),y1(:),'k.');title('Inverse Transform Sampling on rows');
axis equal;axis off;
%% process per column
imcol=cumsum(im,1);
imcol=imcol*size(imcol,2)./repmat(max(imcol,[],1),size(imcol,1),1);
y=1:size(imcol,1);
yimcol_i=zeros(size(imcol));
for i = 1:size(imcol,2)
mask =logical([diff(imcol(:,i))>=0.01;0]);
yimcol_i(:,i) = interp1(imcol(mask,i),y(mask),y);
end
y=1:size(imcol,2);
y=repmat(y,size(imcol,1),1);
y1=y(1:5:end,1:5:end);
yimcol_i1=yimcol_i(1:5:end,1:5:end);
figure('Color','w');plot(y1(:),yimcol_i1(:),'k.');title('Inverse Transform Sampling on cols');
axis equal;axis off;
It has the shortcoming that I can only use this per-row or per-column, but not both. The Inverse Transform Sampling method does not work for multivariate PDFs in general, and I'm fairly sure I wont be able to get it to work in this case.
Is there a simple method to achieve my goal that I haven't seen yet?
I am aware that an algorithm called Voronoi Stippler has been used to create the desired result and I will investigate that, but for the moment I liked the simplicity of Inverse Transform Sampling and would like to know if I can extend that method to match my needs.
It turns out this is fairly simple and can be done by Rejection Sampling.
For the special case where the instrumental distribution is U(0,1) it works like this (if I understood it correctly):
im=imread('http://goo.gl/DDwV3t'); %load random headshot from google
im=imadjust(im,stretchlim(im,[.01,.65]),[]);
im=im2double(rgb2gray(im));
im=im(10:end-5,50:end-5);
im=1-flipud(im);
d = im > .9*rand(size(im));
d=d&(rand(size(d))>.95); %randomly sieve out some more points
[i,j]=ind2sub(size(d),find(d));
figure('Color','w');plot(j,i,'k.');title('Rejection Sampling');
axis equal;axis off;
The sampling is done in one line:
d = im > .9*rand(size(im));
Since I ended up with too many points I randomly sampled the result thus reducing the number of points by approximately the factor 20.
This is pretty much the result I originally desired.