This script demonstrates how to set up a blurring problem. In sequential demos, we will show how the toolbox can solve this problem, even in the presence of noise.
The operator created by this demo is the same created by bbdemo_deblurop, which is used in the remaining demos.
To start it directly, use: playshow bbdemo_deblur1
First we load a test-image so we can see what is going on. This image represents the solution we will seek, represented as a vector X, where the measurement is the transformed image B=A*X (possibly degraded by noise).
To display it properly, we need to compensate for the non-linearity of the display. This is most easily accomplished using BBIMAGEIN, which also shows negative values in blue. We also scale it to make it faster and easier to inspect.
This demo creates an operator, A, which corresponds to a simple model of an imaging system with a stop and an unfocused lens.
dim1=[256,256]; img=double(imread('bird_bw.png')); bb=bbrescale(size(img),dim1,'linear','average'); img=reshape(bb*img(:),dim1); bbimagein(img); axis image off;
Normally, images in Matlab are stored as a matrix or multi-dimensional array. BBTools always assumes that data are represented as a vector, and is usually ignorant about what this vector exactly represents.
To model a "lens stop" we could set a number of pixels to zero. If mask is 1 for pixels we want to keep and 0 elsewhere, the operator could be created as M = bbdiag(mask(:));
However, there is no reason to keep pixels known to be zero. We may therefore use bbindex to create an operator, Mi, that extracts the pixels we want to work further on. The operator M=Mi'*Mi would then act as the diagonal operator above.
To create an image of a circle, we shall use the auxiliary function BBPHANTOM.
mask=bbphantom(dim1,'circle',110); inx=find(mask>.5); Mi=bbindex(inx,prod(dim1)); img1=Mi*img(:); bbimagein(reshape(Mi'*img1,dim1)); axis image off;
A simple, but effective, model of image blurring is a convolution. This mean that each point in the image is splattered out to the shape of a kernel. We will create a "pill-box" kernel (a filled circle), but any shape would be possible here.
The function bbconvn creates an operator that convolves an array (represented as a vector) with a specified kernel. Since we just cut out a circular part, we first need to get the representation back to a rectangular form by multiplying with Mi'.
Also, we can improve performance by cutting out the smallest rectangular part that contains all the data before the convolution. One way to do this is to use bbredim.
Side note: the last step improves performance as it allows the convolution operator to work on a 256-by-256 array. To helps the operator, which is based on the FFT.
K=bbkernel([16,16],'pillbox'); [C,dim2]=bbredim(dim1,[19,19;238,238]); [B,dim3]=bbconvn(dim2,K,'full');
We can now combine all the operators into one: 1) Mask out non-visual part of image (simulate the "stop" of the lens) 2) Cut out the central part 3) Apply blur
To apply this operator, we just have to do a matrix-vector product.
A=B*C*Mi'; b=A*img1; bbimagein(reshape(b,dim3)); axis image off;
As before we have a number of pixels that is known to be zero. Note, however, that the blurring operator increase the number of non-zero pixels in the image.
We do not want to keep around vectors filled with zeros when we get to solving the system. Fortunately we know that the operator contains only non-negative elements (this would not be true of the blurring kernel contained negative values).
A simple way to find the active output-pixels is to multiply by a vector of ones. Rounding-errors means that this must be thresholded to find the true index vector.
f=A*Mi*ones(prod(dim1),1); inx=find(f>=1e-6); Mo=bbindex(inx,size(A,1)); A=Mo*A; bbimagein(reshape(f,dim3)); axis image off;
This completes the blurring operator. Apart from the operator, A, itself, we should also keep the index-operators Mi and Mo to convert between the rectangular form Matlab use for images and the vectors for the problem that only keep non-zero entries.
Without the index operators, Mi and Mo, we would have to deal with a 65536-by-65536 system. Now the system is reduced to 43625-by-38024. This saves quite a bit of effort in the following computations.
The following demos will use this operator, and try to undo it.
plot(0,0); axis([-1.2,1.25,-3,1],'off'); bbshowtree(A,'parent',gca);
