# Documentation

### This is machine translation

Translated by
Mouseover text to see original. Click the button below to return to the English verison of the page.

## Large-Scale Constrained Linear Least-Squares

This example shows how to recover a blurred image by solving a large-scale bound-constrained linear least-squares optimization problem.

### The Problem

We will add motion blur to a photo of Mary Ann and Matthew sitting in Joe's car, then try to restore the original. Our starting image is this black and white image, contained the m x n matrix P. Each element in the matrix represents a pixel's gray intensity between black and white (0 and 1).

```load optdeblur P [m,n] = size(P); mn = m*n; imshow(P); colormap(gray); axis off image; title([int2str(m) ' x ' int2str(n) ' (' int2str(mn) ') pixels'])```

We can simulate the effect of vertical motion blurring by averaging each pixel with the 5 pixels above and below. We construct a sparse matrix D, that will do this with a single matrix multiply.

Create D.

```blur = 5; mindex = 1:mn; nindex = 1:mn; for i = 1:blur mindex=[mindex i+1:mn 1:mn-i]; nindex=[nindex 1:mn-i i+1:mn]; end D = sparse(mindex,nindex,1/(2*blur+1));```

Draw a picture of D.

```cla axis off ij xs = 25; ys = 15; xlim([0,xs+1]); ylim([0,ys+1]); [ix,iy] = meshgrid(1:(xs-1),1:(ys-1)); l = abs(ix-iy)<=5; text(ix(l),iy(l),'x') text(ix(~l),iy(~l),'0') text(xs*ones(ys,1),1:ys,'...'); text(1:xs,ys*ones(xs,1),'...'); title('Blurring Operator D ( x = 1/11)')```

### Apply Blurring Operator

We multiply the image by this operator to create the blurred image. P is the original image, D is the operator, and G is the blurred image.

```G = D*P(:); imshow(reshape(G,m,n));```

### The Blurred Image

Now, let's pretend Joe took this blurred picture G from a moving elevator. Assume we know how fast the elevator is moving, so we know the blurring operator D. How well can we remove the blur and recover the original image P?

The simplest approach is to solve the least squares problem:

```subplot(1,2,1); imshow(reshape(G(:),m,n)); title('G, the blurred image'); subplot(1,2,2); imshow(reshape(P(:),m,n)); title('P, the original image');```

### Regularization Term

In practice the results obtained with this simple approach tend to be noisy. To compensate for this, a regularization term is added:

L is the discrete Laplacian, which relates each pixel to those surrounding it. Since we know we are looking for a gray intensity, we also impose the constraint that the elements of P must fall between 0 and 1.

Create L.

```L = sparse( [1:mn,2:mn,1:mn-1], [1:mn,1:mn-1,2:mn], ... [4*ones(1,mn) -1*ones(1,2*(mn-1))] );```

Draw a picture of L.

```subplot(1,1,1) ; delete(gca); axis ij axis off; xs=11; ys=11; xlim([0,xs+1]); ylim([0,ys+1]); [ix,iy]=meshgrid(1:(xs-1),1:(ys-1)); four=(ix==iy); one=(abs(ix-iy)==1); text(ix(one),iy(one),'-1') text(ix(four),iy(four),'+4') text(ix(~four & ~one),iy(~four & ~one),' 0') text(xs*ones(ys,1),1:ys,'...'); text(1:xs,ys*ones(xs,1),'...'); title('L, The Discrete Laplacian')```

### Set Up Problem to Deblur Image

To obtain the deblurred picture we want to solve for P:

We can simplify this expression by defining A and b:

` A = [D ; 0.02*L]; b = [ G(:) ; zeros(mn,1) ];`

which changes the last equation to:

subject to . Both matrices D and L relate each pixel to a few of its neighbors. This makes A structured and sparse.

```A = [D ; 0.02*L]; b = [ G(:) ; zeros(mn,1) ]; spy(A'); axis equal tight title('Structure of Matrix A''');```

### Solve Optimization Problem

Because A is sparse, we can use a large-scale algorithm to solve this linear least squares optimization problem. We call LSQLIN with A, b, lower bounds, upper bounds, and options.

```options = optimoptions('lsqlin','Algorithm','trust-region-reflective',... 'Display','off'); x = lsqlin(A, b, [], [], [], [], zeros(mn,1), ones(mn,1), [], options); imshow(reshape(x,m,n))```

### Compare Images

Let's compare the blurred and deblurred pictures.

```subplot(1,2,1); imshow(reshape(G,m,n)); title('Blurred'); subplot(1,2,2); imshow(reshape(x,m,n)); title('De-Blurred');```