Computational VisionU. Minn. Psy 5036Linear systems, Convolutions and Optical blurüLinear SystemsLinear system models are important to vision for modeling: optical, retinal sampling, and neural transformations of image data.The basic requirements for a system T to be linear are:T(a+b) = T(a) + T(b)T(λa) = λT(a)For continuous models, the so-called linear superposition operator for an image transformation can be expressed as:üExercises1. As an pencil and paper exercise, show that the above 2D linear superposition operator satisfies the above two criteria.2. Use Mathematica to verify that matrix multiplication on vectors satisfies the two criteria for a linear sytem: superposition andhomogeneity. The matrix WW is the transformation applied to vector inputs aa or bb. aa = Table[a[i],{i,1,3}];bb = Table[b[i],{i,1,3}];WW = Table[W[i,j],{i,1,3},{j,1,3}];üAnswer to 2--or at least one way of doing itSimplify[WW.aa + WW.bb - WW.(aa + bb)]Simplify[WW.(k aa) - k WW.aa]80, 0, 0<80, 0, 0<üConvolutionsOften we can assume that the transformation is space invariant. This is a reasonable assumption if we restrict ourselves to smalloptical patches (so-called isoplanatic patches), uniform retinal sampling, or subgroups of neural systems whose receptive fieldssubtend small regions of the visual field. For a space invariant system, W(x,y; x',y') is equal to W(x-x',y-y') and the linear operation becomes a continuous convolution:from -infinity to +infinity :For computations, we usually model the system operation as a discrete matrix multiplication. Let W be the matrix, r and l arevectors: r = WI2D convolution is used to model: 1) optical blur, 2) discrete retinal sampling by the photoreceptors, and 3) "neural image" transformations by ganglion cell arrays, and cortical cells.Convolution is also used in computer vision as a pre-process for edge-detection, in models of color constancy, and motiondetection. So it is worth spending time to understand it.üContinuous 1-D convolutiona = 3; b = 5;Clear[filter]; (* When you play with building up definitions, you should clear the function to make sure you don't have constraints left hanging from earlier tries*)filter[x_]:= N[(1-Abs[x]/a)] /; x>-a && x <0filter[x_]:= N[(1-x^2/b^2)] /; x < b && x>=0filter[x_]:= 0.0 /; x<=-a || x >= bPlot@filter@xD, 8x, -8, 8<, PlotRange Ø 8-1, 1<D-55-1.0-0.50.51.0/; means "such that". Mathematica uses notation borrowed from the C programming language for logical operations such as: OR => || and AND => &&edge[x_] := N[If[x<4,1/2,1.5]];2 Convolutions_Tutorial2.nbPlot@edge@xD, 8x, -10, 10<D-10-55100.81.01.21.4A 1-D convolution operation is written:It is useful to visualize this as the area underneath the curve formed by the product of edge and the left-right reversed filter as itslides along the x-axis.Animate@Plot@edge@xD filter@z - xD, 8x, -30, 30<, PlotRange Ø 88-15, 15<, 8-1, 2<<D, 8z, -10, 10, 0.5<Dz-15-10-551015-1.0-0.50.51.01.52.0We can find the area under each of the above curves by using Mathematica's built-in numerical integration capability, NInte-grate[]. It helps if we specify the range of necessary integration as precisely as possible. Inspection of the above curves showsthat {z-b,z+a} is an appropriate range for x.Convolutions_Tutorial2.nb 3r1 =Table[{z,NIntegrate[edge[x] filter[z - x],{x,z-b,z+a}]}, {z,-20,20,1}];NIntegrate::slwcon :Numerical integration converging too slowly; suspect one of the following: singularity, valueof the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. àNIntegrate::slwcon :Numerical integration converging too slowly; suspect one of the following: singularity, valueof the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. àcontinousconvg = ListPlot@r1, PlotRange Ø 8-1, 8<, PlotStyle Ø 8RGBColor@1, 0, 0D<D-20-1010202468üDiscrete 1-D convolutionIn this section, we will pass the 1-D "edge" through the filter again, but this time we will set it up as a discrete matrix operation.We'll apply a convolution matrix to a vector that represents a sharp edge, l=edge, at x = size/2. size = 30; edge@x_D := NBIfBx < 4,12, 1.5`FF;edgelist = TableB8x, edge@xD<, :x, -size2,size2- 1>F; edgevector = Transpose@edgelistDP2T;edgelistg = ListPlot@edgelist, PlotRange Ø 8-1, 8<, PlotStyle Ø 8RGBColor@0, 0.5`, 0D<D-15-10-55102468Because the integration is finite, we have to make some decision about what to do at the boundaries. One choice is to assume thatw and l are zero beyond the boundaries. Another is to assume that they are periodic with a shared common period, and we are justrepresenting one period of the signal. Now we construct a size x size matrix that will perform a convolution operation.4 Convolutions_Tutorial2.nbBecause the integration is finite, we have to make some decision about what to do at the boundaries. One choice is to assume thatw and l are zero beyond the boundaries. Another is to assume that they are periodic with a shared common period, and we are justrepresenting one period of the signal. Now we construct a size x size matrix that will perform a convolution operation.w = Table[filter[i-j], {i,size},{j,size}];A convenient way to view the matrix is to use ArrayPlot[], where you can see how each row is a shifted version of the precedingrow. ArrayPlot@w, ImageSize Ø TinyDrvect = w.edgevector; r = TableB:i -size2, rvectPi + 1T>, 8i, 0, size - 1<F; matrixconvolg =ListPlot@r, Joined Ø True, PlotRange Ø 8-1, 8<, PlotStyle Ø 8RGBColor@0, 0, 0.5`D<D;Show@8matrixconvolg, edgelistg, continousconvg<D-20-1010202468The convolution matrix "blurs" out the edge and is in reasonable agreement with the approximation to the continous convolution.Note that the left and right boundary edges also got "blurred" by the matrix convolution. üExercise: Blur the above edge with a line-spread function (LSF) that models diffraction-limited opticsüIntroduction to diffraction limited blurPin-hole optics are limited by diffraction. The image of a point though a circular aperture is an Airy pattern. We have to take careto define the function at x and y =0.(*Cell inactive*)Airy2D[x_,y_] := If[x==0 && y==0,1,(2 BesselJ[1,Pi Sqrt[x^2+y^2]]/(Pi Sqrt[x^2+y^2]))^2]Plot[Airy2D[x,0],{x,-3,3}];For a single slit aperture, the "line-spread-function" for a diffraction limited system is a "sync" function squared: Convolutions_Tutorial2.nb 5For a single slit aperture, the