Random Image Segmentation

I’ve been playing around today with an image decomposition method. Given some observed image \bf y, one seeks an image \bf x so as to minimize

\sum_i (x_i - y_i)^2 + \lambda \sum_{(i,j)} |x_i-x_j|

Where the first sum is over all pixels i, and the second sum is over all neighboring pairs of pixels (i,j) (For simplicity here, I imagine that images are just long vectors, and the grid structure is only present in what pairs are considered neighbors). This minimization has the interesting property that, at a solution, neighboring pixels in \bf x often have the same value. Thus, one can see the optimization as segmenting the image into piecewise constant regions.

I tried using the (excellent) CVX toolbox to do the optimization. Though this worked, it was quite slow. I then derived a little iteratively reweighted least-squares algorithm that is much faster.

(Update Jan 9, 2009: Code for doing the optimization is here.)

Anyway, the standard use of the decomposition is for “denoising”. Here is an example from an image in the Street Scenes database. At top is the original image \bf y, followed by the result \bf x. At the bottom, I use a random colors for each segment to give a better idea of the discontinuities.

Always, intensities are in the 0-1 interval, and only immediate horizontal and vertical pairs are considered neighbors.

\lambda=1

\lambda=1

However, more interesting (to me) images come from segmenting random images.

Uniform noise, \lambda=0.1

Uniform noise, \lambda=0.2

Uniform noise, \lambda=0.3

Uniform noise, \lambda=0.4

Uniform noise, \lambda=0.5

Uniform noise, \lambda=0.6

Gaussian noise, variance \sigma^2=.1, \lambda=0.3

Gaussian noise, variance \sigma^2=.1, \lambda=0.4

Gaussian noise, variance \sigma^2=.1, \lambda=0.5

Advertisements

Markov’s Inequality

While looking at Tao’s post on the law of large numbers, I claimed to Alap that Markov’s inequality was “obvious”. Asked to back that up, though, I couldn’t remember why it was obvious. Take some distribution p(x). For simplicity, assume this distribution is over the non-negative reals. Then, Markov’s inequality states, for \lambda > 0

P[x\geq\lambda] \leq \frac{E[x]}{\lambda}.

Though I find it more intuitive in the form

\lambda P[x \geq \lambda] \leq E[x].

Now, a proof is pretty trivial, basically by writing down the definitions of both sides.

E[x] = \int_{x\geq 0} x p(x) dx

\lambda P[x \geq \lambda] = \int_{x\geq\lambda} \lambda p(x) dx \leq \int_{x\geq\lambda} x p(x) dx

\leq \int_{x\geq 0} x p(x) dx

However, I think more intuition can be conveyed with pictures.

Matlab Code:

function markov_inequality

fsize = 16;
lambda = 5;

x = 0:.01:15;
p = .5*exp(-(x-3).^2) + exp(-(x-7).^2/4);
p=p/sum(p)*100;

figure(1)
plot(x,p,'k-','LineWidth',2)
xlim([min(x) max(x)])
ylim([0 max(x.*p)+.25])
legend('p(x)');
set(gca,'FontSize',fsize)

figure(2)
plot(x,p.*x,'k-','LineWidth',2);
xlim([min(x) max(x)])
ylim([0 max(x.*p)+.25])
legend('x p(x)');
set(gca,'FontSize',fsize)

figure(3)
plot(x,p.*x,'k-',x,p.*lambda,'r-','LineWidth',2);
xlim([min(x) max(x)])
ylim([0 max(x.*p)+.25])
hold on;
legend('x p(x)', '\lambda p(x)');
title('\lambda=5');
set(gca,'FontSize',fsize)

figure(4)
area(x,x.*p,'FaceColor',[.7 .7 .7],'LineWidth',2)
hold on;
good = find(x>=lambda);
area(x(good),p(good)*lambda,'FaceColor',[.8 0 0],'LineWidth',2)
hold off
legend('E[x]','\lambda P[x \geq \lambda]')
xlim([min(x) max(x)])
ylim([0 max(x.*p)+.25])
set(gca,'FontSize',fsize)

end

More Implicit Functions

(This time, including derivatives.) Closed form solutions are possible for some of these, I know.






…which of course has a constant derivative.

Matlab code:

function implicit_funcs
options = optimset('TolX',1e-15,'MaxFunEvals',1e5,'MaxIter',1e5);
fsize = 16;

e = 1e-5;

f1 = @(x) 5*x;
f2 = @(x) x;

xAR = -3:.05:3;
for n=1:length(xAR)
    x = xAR(n);
    myfun = @(y) 5*abs(y-f1(x))^2 + abs(y-f2(x))^2;
    y_star(n) = fminsearch(myfun,0,options);
    
    % approximate derivative with one sided finite diff
    x = xAR(n)+e;
    myfun = @(y) 5*abs(y-f1(x))^2 + abs(y-f2(x))^2;
    y_star2(n) = fminsearch(myfun,y_star(n),options);
end
figure(1)
plot(xAR,f1(xAR),'k-',xAR,f2(xAR),'k--',xAR,y_star,'r-','LineWidth',2)
legend('f1','f2','y*')
title('y*(x) = arg min_y = |y-f1(x)|^{2}+ (y-f2(x))^2','FontSize',fsize)
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)
set(gca,'FontSize',fsize)

figure(2), plot(xAR,(y_star2-y_star)/e,'r-','LineWidth',2);
title('y*(x) = arg min_y = |y-f1(x)|^{2}+ (y-f2(x))^2','FontSize',fsize)
legend('dy*(x)/dx')
set(gca,'FontSize',fsize)
xlabel('x','FontSize',fsize)

end

Implicit Functions


Matlab code:

function implicit_funcs

fsize = 16;

f1 = @(x) x;
f2 = @(x) 1.5*x;

xAR = -2.5:.05:2.5;
for n=1:length(xAR)
    x = xAR(n);
    myfun = @(y) abs(y-f1(x)) + (y-f2(x))^2;
    y_star(n) = fminsearch(myfun,0);
end
plot(xAR,f1(xAR),'k-',xAR,f2(xAR),'k--',xAR,y_star,'r.','LineWidth',2)
legend('f1','f2','y*')
title('y*(x) = arg min_y = |y-f1(x)|+ (y-f2(x))^2','FontSize',fsize)
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)
set(gca,'FontSize',fsize)

clear
figure
fsize = 16;

I=3;
f{1} = @(x) x;
f{2} = @(x) 1.5*x+1;
f{3} = @(x) .5*x+1;
function o=myfun2(y)
    o = 0;
    for i=1:I
        o = o + abs(y-f{i}(x));
    end
end
xAR = -5:.1:5;
for n=1:length(xAR)
    x = xAR(n);
    y_star(n) = fminsearch(@myfun2,0);
end
plot(xAR,f{1}(xAR),'k-',xAR,f{2}(xAR),'k-',xAR,f{3}(xAR),'k-',xAR,y_star,'r.','LineWidth',2)
legend('f1','f2','f3','y*')
title('y*(x) = arg min_y = \Sigma_i |y-f_i(x)|','FontSize',fsize)
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)
set(gca,'FontSize',fsize)

clear
figure
fsize = 16;

I=3;
f{1} = @(x) x;
f{2} = @(x) 1.5*x+1;
f{3} = @(x) .5*x+1;
function o=myfun3(y)
    o = 0;
    for i=1:I
        o = o + (y-f{i}(x))^2;
    end
end
xAR = -5:.1:5;
for n=1:length(xAR)
    x = xAR(n);
    y_star(n) = fminsearch(@myfun3,0);
end
plot(xAR,f{1}(xAR),'k-',xAR,f{2}(xAR),'k-',xAR,f{3}(xAR),'k-',xAR,y_star,'r.','LineWidth',2)
legend('f1','f2','f3','y*')
title('y*(x) = arg min_y = \Sigma_i (y-f_i(x))^2','FontSize',fsize)
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)
set(gca,'FontSize',fsize)

end