Convolution and Fourier Transformation PDF

Convolution and Fourier Transformation Segmentation & Visualization Ferdi van der Heijden  University of Twente - RAM  10/28/2021 Contents 1 INTRODUCTION.............................................................................................................................

Convolution and Fourier Transformation Segmentation & Visualization Ferdi van der Heijden  University of Twente - RAM  10/28/2021 Contents 1 INTRODUCTION........................................................................................................................................... 1 2 POINT SPREAD FUNCTIONS.......................................................................................................................... 1 2.1 THE POINT SPREAD FUNCTION IN OPTICAL SYSTEMS................................................................................ 2 2.2 THE POINT SPREAD FUNCTION IN X-RAY SYSTEMS................................................................................... 3 2.3 OTHER POINT SPREAD FUNCTIONS........................................................................................................ 4 2.4 PROPERTIES OF POINT SPREAD FUNCTIONS............................................................................................ 4 3 2D CONVOLUTION...................................................................................................................................... 5 3.1 THE SIFTING INTEGRAL....................................................................................................................... 5 3.2 THE CONVOLUTION INTEGRAL............................................................................................................... 7 3.3 EDGE AND LINE SPREAD FUNCTIONS..................................................................................................... 9 3.3.1 THE LINE SPREAD FUNCTION........................................................................................................... 9 3.3.2 THE EDGE SPREAD FUNCTION....................................................................................................... 10 3.4 DISCRETE CONVOLUTION................................................................................................................. 11 3.4.1 IMPLEMENTATION DETAILS........................................................................................................... 11 3.4.2 PROPERTIES OF CONVOLUTION..................................................................................................... 13 3.5 APPLICATIONS................................................................................................................................ 14 3.5.1 IMAGE ENHANCEMENT AND NOISE REDUCTION................................................................................ 14 3.5.2 IMAGE FEATURE EXTRACTION........................................................................................................ 15 4 2D FOURIER TRANSFORMATION................................................................................................................ 17 4.1 THE FOURIER INTEGRAL................................................................................................................... 17 4.1.1 COMPLEX HARMONICS................................................................................................................ 17 4.1.2 PROPERTIES OF COMPLEX HARMONICS.......................................................................................... 19 4.1.3 THE FOURIER TRANSFORMATION.................................................................................................. 20 4.1.4 EXAMPLES................................................................................................................................. 23 4.2 OPTICAL TRANSFER FUNCTIONS........................................................................................................ 24 4.3 DISCRETE FOURIER TRANSFORM....................................................................................................... 26 4.3.1 DISCRETIZATION OF THE FOURIER TRANSFORMATION....................................................................... 26 4.3.2 PROPERTIES.............................................................................................................................. 27 4.4 IMAGE FILTERING WITH THE DFT....................................................................................................... 29 4.4.1 FROM PSF TO OTF.................................................................................................................... 30 4.4.2 BOUNDARY EFFECTS................................................................................................................... 30 4.5 APPLICATIONS................................................................................................................................ 31 4.5.1 IMAGE TEXTURE EXTRACTION........................................................................................................ 32 4.5.2 IMAGE RESTORATION................................................................................................................... 32 APPENDIX A PROPERTIES OF 2D FOURIER TRANSFORMS................................................................................... 36  1 Convolution and analysis and synthesis. This lecture note is an introduction to the topic. Being the Fourier pillars of linear image filtering, we introduce convolution and Fourier analysis. Convolution is first presented from a Transformation physical point-of-view. This starts with an introduction of the point spread function leading to convolution in the physical, continuous imaging domain. Next, the transition to the Segmentation & Visualization discrete, digital domain is made. The result is discrete convolution, suitable for implementation in MATLAB and Python. 1 Introduction Fourier analysis will also be introduced in the continuous One of the more important operations in digital image domain. After this, the switch to the discrete domain will be processing is image filtering. This operation has many made. The relationship between Fourier transformation and applications: image reconstruction, image enhancement, convolution is the next topic. This lecture note concludes image feature extraction, image restoration, etc. with a demonstration of how Fourier analysis in the frequency domain can be helpful for the implementation of Image filtering is a neighbourhood-to-pixel operation. This is convolution in contrast with grey level transformations which are pixel- to-pixel operations: an output pixel depends only on the 2 Point spread functions corresponding input pixel. In image filtering, an output pixel depends on a neighbourhood of the corresponding input Figure 2 is an illustration of the camera obscura. This pixel. See Figure 1. This also works the other way round: an device is regarded as a forerunner in photography. It is not input pixel has its influence in the neighbours of the much more than a closed box with on one side a tiny hole, corresponding output pixel. called the pinhole. Image filtering is either linear or nonlinear. The nonlinear Under the assumption that light travels along straight lines, operations have in principle more potentials. Nevertheless, any 3D light emitting point in front of the box projects its linear image filtering is important because of its light along a ray through the pinhole and is seen as a small mathematical tractability, which enables systematic spot on the opposite wall in the box. This is schematically shown in Figure 3. If X = ( X , Y , Z ) is the coordinates of a 3D point in space, with a coordinate system as indicated in the figure, then this point is projected on the 2D image plane at a position: Xd Yd x0 = − y0 = − (1) Z Z Figure 1. Neighbourhood operation. Top: an output pixel gn,m depends on a neighbourhood of the corresponding input pixel fn,m. Bottom: at the same time, an input pixel Figure 2. Engraving of a portable camera obscura. influences the output image in a neighbourhood. Athanasius Kircher: Ars Magna Lucis Et Umbrae (1645)  2 Figure 4. Aperture function The important phenomenon to observe is that, apart from Figure 3. Projection of a point light source in a camera the factor, the image of the point source only shifts: it does obscura. not change shape. The image of the point source is called The direction of Z is defined orthogonal to the image. d is the point spread function, abbreviated to PSF. the distance between the pinhole and the image plane. It is called the focal distance. The projection that is expressed 2.1 The point spread function in optical systems in eq (1) is called the perspective projection. Image formation in optical cameras follows the model of a pinhole. See Figure 5. The aperture is replaced by a lens A 3D point source is not imaged as a point. This is because system. This reduces the spatial extent of the point spread the pinhole must have a finite size. See Figure 4. If it would function, but still it cannot prevent that a point source is be infinitely small, no light energy can enter the box, and imaged as spot rather than a point. Figure 5 shows the point would be invisible. measured PSFs of six lenses. For comparison of the size of A finite size of the pinhole implies that the 3D point source the PSFs, the white square has a size of 8.5 µ m which is is imaged by a spot rather than a point. The aperture typical for the pixel size of 12 MP cameras. function of the system describes the shape of the pinhole, The PSF of images 1 to 5 are suitable for moderate and is mathematically defined by: cameras. These PSFs all exceed the size of the pixel. The 1 if ( x, y ) is inside the pinhole PSF shown in image 6 is from an excellent image quality, h ( x, y ) =  (2) 0 if ( x, y ) is outside the pinhole When projected to the image plane, the aperture is magnified by a factor = s ( Z + d ) Z. If the 3D point source is lying far away from the pinhole, Z >> d , which we assume, then s ≈ 1. Suppose that the 3D point is lying on the Z-axis, so that its position is X = (0, 0, Z ) , then the image of the point source is: I0 f ( x, y ) = h ( x, y ) (3) Z2 I 0 is the intensity of the point source. The factor 1 Z 2 accounts for the increase of spread of radiant energy as Z increases. If the point source shifts to an arbitrary 3D position with coordinates X = ( X , Y , Z ) , its image shifts to ( x0 , y0 ) , given in eq (1). The image of the 3D point becomes: I0 f ( x, y ) cos3ϕ = h ( x − x0 , y − y0 ) (4) Z2 The factor cos3ϕ accounts for the apparent foreshortening Figure 5. The point spread functions (PSF) of optical lenses. From: H.H. Nasse “How to Read MTF Curves”, of the aperture when illuminated from a sideward position. Carl Zeiss AG (2008).  3 and its size is smaller than the pixel size. For the 12 MP from the sign, the projection is the same as in the pinhole camera, this lens would be over specified. model expressed in eq (1). 2.2 The point spread function in X-ray systems Because of the extreme concentration of heat energy in a single point it is difficult to create a true X-ray point source. Figure 6 shows an X-ray system consisting of an X-ray In the real physical world, X-rays originate from a finite flat source and an X-ray detector, in this case an image area. Not every part in this area contributes equally to the intensifier. The geometry of the image formation is also total X-ray energy. The total radiant exitance is distributed shown in the figure. over the area according to a distribution function E ( x, y ) : An ideal X-ray source is a point source. Objects between the focal spot distribution. For pure technical reasons, the source and image plane attenuate the X-rays so that a distribution has two elongated ridges. Each ridge has a size shadow image is created. If the object is just a point, the of the order of magnitude of 0.3 × 1 mm. Next to each other, shadow would also be a point in the image plane. Assuming the total size of the two ridges is roughly around 1× 1 mm. that a coordinate system has its origin at the point source, If a point-shaped object, that absorbs the X-ray radiation, is and the z-direction is orthogonal to the image plane, a placed at a position X = ( X , Y , Z ) , then this creates an point-shaped object at position X = ( X , Y , Z ) will be image of the focal spot distribution at the image plane. This imaged at: image is proportional with (Figure 7): Xd Yd =x0 = y0 (5) Z Z Z f ( x, y ) = E ( s ( x0 − x), s ( x0 − x) ) with s = (6) Z −d d is the distance between source and image plane. Apart Thus, the PSF of the system is: h( x, y, Z ) = E (− sx, − sy ) (7) In Figure 7, the focal spot plane is parallel with the image plane. As a result, the PSF in eq (7) does not depend on ( x0 , y0 ). With fixed Z , the PSF is space invariant. In real life, the situation is more complicated. The focal spot plane is not parallel with the image plane; it is tilted. This is done to enlarge the area at which heat can be dissipated. The consequence is that, seen from the image plane, the focal spot area is foreshortened, so that the PSF looks smaller. But the price to pay is that the PSF now becomes Figure 6. X-ray fluoroscopy system. Top: Philips C-arm. Bottom: perspective projection in such system. Figure 7. The point spread function in an X-ray system.  4 Figure 8. Measured point spread functions of an X-ray Figure 9. The point spread function of a defocused camera. system. Their shapes are space-variant. Top: image plane a is in focus; plane b is out of focus. From: van der Heijden, F, 1981 Bottom: the point spread function at plane b is a disk. space variant. The shape changes with varying ( x0 , y0 ).  x2 + y 2  h= ( x, y ) A exp  −  (10)  2σ 2  Figure 8 shows the result of an experiment. A leaden plate was constructed with tiny holes in it at a regular spacing. Point spread functions are also used to model motion blur. The holes are small compared with the size of the focal Figure 12 is a long exposure photo of a starry night. On the spot. Each hole acts like a camera obscura, thus showing scale of this picture, the stars are point sources. The an image of the focal spot distribution. Figure 8 shows motion of the stars causes them to be trajectories in the clearly the space variant shape of the distribution. image plane and thus to form streaks. These streaks are the PSFs. As the motion is not uniform, the functions are 2.3 Other point spread functions space variant. In many imaging systems, the PSF plays a major role in the image quality. Examples are photographic paper and film, 2.4 Properties of point spread functions display units and projectors, acoustic imaging, thermal Rotational symmetry imaging, X-detectors such as image intensifiers, etc. A point spread function h( x, y ) is fully rotational symmetric A daily life example is the defocused camera. See Figure 9. if a rotation of the coordinate system at any angle does not The PSF is approximately spatial invariant. It has the shape change the PSF. Let this rotation be described by of a disk. If R is the radius of the disk, the mathematical = ξ x cos α + y sin α and = η x sin α − y cos α , then the PSF expression is:  A if x 2 + y 2 ≤ R 2 h ( x, y ) =  (8) 0 elsewhere The 2D rectangular function is another PSF, defined by:  A if − 12 X < x ≤ 12 X and − 12Y < y ≤ 12Y h ( x, y ) =  (9) Figure 10. Rectangular point spread function. 0 elsewhere See Figure 10. This function is used to model the behaviour of image sensors, where each pixel catches photons that fall on a rectangular area. Yet another PSF is the Gaussian, shown in Figure 11: Figure 11. Gaussian point spread function.  5 The rectangular function in eq (9) and the Gauss function in eq (10) are examples of separable functions. The volume of a point spread function The graphs in figures 9, 10 and 11 show that the PSF can be regarded as 2D surfaces in a 3D space. The volume of the space below the surface and the base plane is a volume: +∞ +∞ V= ∫ ∫ h( x, y)dxy −∞ −∞ (13) This volume is one of the characteristics of a PSF. It is easy Figure 12. Point spread functions associated with motion to calculate this volume for the three given cases: blur. Photo: F. Rødland Disk function: V = Aπ R 2 is rotational symmetric if for any rotation angle α , we have Rectangular function: V = AXY h(ξ ,η ) = h( x, y ) Gauss function: V = 2π Aσ 2 The function is rotationally symmetric if and only if it is only Homework 2_1: a function of the radius=r x2 + y 2 : Make a sketch of the 1D function h(r ) of the disk function and the Gauss function. h( x,= y ) h(r , 0) with=r x2 + y 2 (11) Give the mathematical expression for h(r ) of the disk function and the Gauss function. The 2D function can then also be represented by a 1D Make a sketch of the two 1D functions h1 ( x) and function: h(r ) = h(r , 0). h2 ( y ) of the rectangular function and the Gauss function. Give mathematical expressions for these The rotational symmetry of a PSF makes it easier to analyse functions its behaviour. For instance, to measure the PSF of an imaging device, it suffices to measure the so-called line 3 2D Convolution spread function in just one direction. See Section 3.3.1. Measuring the line spread function is much easier than The previous section considered the image of a scene that measuring the PSF. From the observed line spread consisted only of a single point source, or a single point function, the PSF can be derived. object. The current section addresses the question how the image would look like for real scenes. The disk function in eq (8), and the Gauss function in eq (10), are both rotationally invariant. 3.1 The sifting integral Separability The proposition of the current section is that: A PSF is separable in x and y if and only if the function "Any image can be regarded as the sum of an can be written as the product of two 1D functions, one infinite number of points". which only depends on x , and one that only depends on y : The sifting integral is the mathematical basis of this statement. h( x, y ) = h1 ( x)h2 ( y ) (12) The Dirac function If a function is separable, it facilitates the analysis. More It starts with the introduction of a 2D Dirac function δ ( x, y ) important is, that if a PSF is separable, the convolution (see. This function is defined such that it holds the following two Section 3) can be implemented much more efficient. That properties: is, in digital image filtering, convolution with a separable PSF requires much less computations than a non-separable one.  6 Figure 13. The Dirac function is the limiting case of a Gauss function as σ → 0. = a) δ ( x, y ) 0 if ( x, y ) ≠ (0, 0) +∞ +∞ (14) Figure 14. An image consisting of two points. b) ∫ ∫ δ ( x, y)dxdy = 1 −∞ −∞ This is demonstrated in Figure 14. The total radiant energy One of the functions that satisfies these requirements is is now a0 + a1. the limiting case of a Gauss function: The Comb function  1  x2 + y 2  Consider an orthogonal grid with a spacing of ∆. Grid =δ ( x, y ) lim  exp  −   (15) σ → 0 2πσ  2  2σ 2  points in x -direction are enumerated by integers n , so that their coordinates are xn = n∆. Grid points in y - This is illustrated in Figure 13. The more σ approaches direction are enumerated by m. Thus ym= m∆. zero, the more the Gauss function satisfies the requirements of a Dirac function. In the limiting case, the Suppose that, on each grid point, a Dirac function is function represents a single point. To guarantee that the situated with intensity f n , m ∆ 2. This defines a so-called function still carries energy, the amplitude approaches comb function: infinity. = f ( x, y ) ∑∑ f n m n,m δ ( x − n ∆ , y − m∆ ) ∆ 2 (18) Suppose that we would have a camera with an ideal lens. If the scene would consist of a single point source, then this Figure 15a shows an example of a comb function defined lens projects this source to a point on the image plane of on a 5x5 grid. Figure 15b shows the same image plane as the image sensor. If the position of this point is ( x0 , y0 ). in Figure 15a. However, the spacing ∆ of the grid is 4x then the mathematical representation of this point in the smaller, so that the area is now populated with a 20x20 image is: grid. Figure 16 shows the comb function in the image plane on an 80x80 grid. f ( x, y ) = a0δ ( x − x0 , y − y0 ) (16) In eq (18), the intensities of the points are f n , m ∆ 2. The The variable a0 is the intensity of the Dirac function, but it reason of embedding the factor ∆ 2 is because the density also quantifies the total radiant energy ∫ ∫ +∞ −∞ f ( x, y ) dxdy on (= number of points per unit area) is proportional to 1 ∆ 2. the image plane. By embedding the factor ∆ 2 , the average radiant energy per unit area does not depend on ∆ As said before, a0δ ( x − x0 , y − y0 ) is the ideal image of a point source, but for brevity’s sake we will call this Dirac When increasing the density of points in the image plane, function itself just a “point”. If we have two of these points the result approaches more and more a real image. In the in an image, the representation is: limiting case, as ∆ approaches zero, the number of points per unit area goes to infinity, while ∆ 2 goes to zero. The f ( x, y )= a0δ ( x − x0 , y − y0 ) + a1δ ( x − x1 , y − y1 ) (17) discrete number of points becomes a continuum.  7 a) b) Figure 16. The comb function defined on an 80x80 grid. Figure 15.Comb functions. a) defined on 5x5 grid, b) defined on a 20x20 grid.   =f ( x, y ) lim  ∑∑ f n , mδ ( x − n∆, y − m∆)∆ 2  (19) ∆→ 0  n m  Figure 17. The continuous image. +∞ +∞ The double summation is a Riemann sum, and in the limiting case it becomes a Riemann integral. By substitution = f ( x, y ) ∫∫ −∞ −∞ f ( x, y )δ ( x − ξ , y − η )d ξ dη (22) of d ξ dη = ∆ 2 , and ξ = n∆ , η= m∆ , and f (n∆, m∆) =f n , m , we symbolically write: 3. Using the fact that δ ( x − ξ , y − η ) is zero everywhere, except at ξ = x and η = y , we can ignore all values of ∞ ∞ f ( x, y ) , except the one with ξ = x and η = y. That is, = f ( x, y ) ∫ ∫ ξ = −∞ η = −∞ f (ξ ,η )δ ( x − ξ , y − η )dη d ξ (20) except for f (ξ ,η ). This brings us to eq (20). QED. This is the sifting integral. It has this name as it acts like a By substitution of ξ= x − ξ and η= y − η , so that sieve. It takes the input image f (ξ ,η ) and selects (sifts) d ξ = −d ξ and dη = −dη , it is easy to prove that the sifting the element from f (ξ ,η ) at a position that is indicated by integral also can be written as: = the position of the Dirac function, i.e., at ξ x=, η y. ∞ ∞ The validity of the sifting integral can also be made f= ( x, y ) ∫ ∫ ξ = −∞ η = −∞ f ( x − ξ , y − η )δ (ξ ,η )dη d ξ (23) plausible in a more mathematical way. Consider the image f ( x, y ) : 3.2 The convolution integral 1. Step 1 is to multiply this function with Now that it is known that an image can be regarded as an ∫ ∫ δ ( x − ξ , y − η )d ξ dη. This is allowed since this infinite collection of points, it is easy to describe what integral is by definition one: happens if an imaging system has a PSF of h( x, y ). We just have to replace the Dirac functions by the PSF. This is +∞ +∞ f ( x, y ) f ( x, y ) ∫ shown in Figure 18 for the case of an input image = ∫ δ ( x − ξ , y − η )dξ dη −∞ −∞ (21) consisting of two points. The PSF is a Gaussian as given in eq (10). In general, the response on the two points is: 2. Since the integral integrates over ξ ,η , and not over x, y , the function f ( x, y ) is constant with respect to g ( x, y )= a0 h( x − x0 , y − y0 ) + a1h( x − x1 , y − y1 ) (24) ξ ,η. Step 2 is to bring this constant within the integral: The addition in this equation is allowed because of the additive nature of light. In Figure 18 the two individual  8 a) b) Figure 18. The response of a system with a Gaussian point spread function on an input image with two points. Figure 19. The response of a system with a Gaussian PSF. responses of the two points overlap. The response of the a) Input is a comb function with a 20x20 grid. b) Input is a continuous image. sum is the sum of the two individual responses. For the comb function, the response on a system with a there are two operands: in our case two input images. PSF h( x, y ) is (see eq (18)): There is only one output. Symbolically, the operation is written as: = g ( x, y ) ∑∑ f n,m h ( x − n ∆ , y − m∆ ) ∆ 2 (25) n m g ( x, y ) = f ( x, y ) * h ( x, y ) (28) Figure 19 shows the output of a system with the Gaussian The operator * denotes convolution and should not be PSF applied to the 20x20 comb function displayed in Figure confused with multiplication as in a computer language. As 15b. It can be seen that, due to the overlapping Gaussian eq (27) shows, the operator is commutative: function, the image is smoother. Actually, for human perception the filtered version is more informative than the f ( x, y ) * h ( x, y ) = h ( x, y ) * f ( x, y ) (29) original spiky one. The used notation underlines the fact that a PSF can also When the imaging system is applied to a continuous image, be regarded as an image. the sifting integral comes into play. Being the Riemann sum of Dirac functions, we just need to replace these Dirac In functional block diagrams, operators are drawn as functions in eq (20) by the PSF. This brings the 2D rectangles, whereas the operands are arrows. Figure 20 convolution integral: provides examples. ∞ ∞ Homework 2_2: = g ( x, y ) ξ ∫ η∫ = −∞ = −∞ f (ξ ,η )h( x − ξ , y − η )dη d ξ (26) Suppose the input image is just a constant: f ( x, y ) = C. Give an expression of the response if the PSF is a disk The result of the convolution operator with a Gaussian PSF, function as in eq (8). applied to the continuous image in Figure 17, is given in Give an expression of the response if the PSF is a Figure 19b. The output image is a smoothed version in rectangular function as in eq (9). which possible noise has been suppressed, but at the same Give an expression of the response if the PSF is a time all image details are fully lost. This behaviour of the Gaussian function as in eq (10). convolution depends heavily on the choice of the PSF. Hint: substitute f ( x, y ) = C in eq (26) or (27). Choose the Other PSF can have an opposite effect. equation which suits you best. Then, work out the integral. Since there are two versions of the sifting integral, see eq (23), there are also two versions of the convolution integral: ∞ ∞ g= ( x, y ) ∫ ∫ ξ = −∞ η = −∞ f ( x − ξ , y − η )h(ξ ,η )dη d ξ (27) Figure 20. Functional block diagram of Convolution is a dyadic operation. “Dyadic” means that the convolution operator.  9 Figure 22. Top view of the disk function. ∞ ∞ = lsf ( x) ∫ ∫ ξ = −∞ η = −∞ δ ( x − ξ )h(ξ ,η )dη d ξ ∞ ∞ ∞ (31) =∫ ∫ δ ( x − ξ )h(ξ ,η )d ξ dη =∫ h( x,η )dη η = −∞ ξ = −∞ η = −∞ In the first step, the order of integration is swapped. In the Figure 21. The line spread function of disk function. second step, the 1D version of the sifting integral has been a) Input image, b) Output image. applied. 3.3 Edge and line spread functions Equation (31) is the relation between point and line spread The PSF and the convolution integral provide good models function. If the PSF is rotational symmetric, then this to describe behaviour of physical imaging systems. With solution is valid for any direction of the wire. If the PSF is knowledge of the PSF, the response of a system to input not rotational symmetric, we have to repeat the procedure structures other than points can be given. Two examples for other directions. are the line spread function and the edge spread function. Example: the line spread function of a disk 3.3.1 The line spread function Consider the PSF given in eq (8) and shown in Figure 9. What will be the LSF? To answer this question, we have to Measuring the PSF of physical imaging system is not easy. substitute h( x,η ) in eq (31), and work this out. The The PSF is often just a tiny spot, which might be problem here is that the disk function is piecewise overwhelmed by noise and hardly be visible. A better continuous: at the circle with radius R it jumps from zero method is to measure the LSF (line spread function). From to A. This is shown in Figure 22. Outside the circle, the a set of LSFs, the PSF can be deferred. For that, the disk function is zero. Inside, it is A. relationship between point and line spread function is needed. To solve ∫ h( x,η )dη , we need to distinguish different intervals of x : Measuring the LSF is done by imaging a thin structure. In X- If x < − R of x > + R , we are always outside the circle. ray imaging, a thin metal wire, which is spanned parallel In that case, h( x,η ) ≡ 0 , and thus ∫ h( x,η )dη is zero. with the image plane, will do. Such a thin object is If − R < x < + R , the function h( x,η ) is piecewise mathematically represented by a Dirac function. Suppose continuous: that the wire crosses the origin and the direction of the wire is along the y-direction, then the model of the wire is:  A if η 2 < η < η1 h( x,η ) =  with η1.2 = ± R2 − x2 (32)  ( x, y ) = δ ( x ) (30) 0 elsewhere See Figure 21a. We can substitute this in eq (26) or eq (27) In this particular case, the integral is:. The latter seems to provide the easiest form here. +∞ η2 η1 +∞ Substitution yields: ∫ −∞ h( x,η )dη = ∫ −∞ h( x,η )dη + ∫ h( x,η )dη + η2 ∫ h( x,η )dη η1  10 η2 η1 +∞ = ∫ 0dη + ∫ Adη + −∞ η2 ∫ 0dη η1 (33) = (η1 − η 2 ) A The line spread function becomes: 0 if x < − R  = lsf ( x) 2 A R 2 − x 2 if − R < x < + R (34) 0 if x > + R  A graphical representation in Figure 21 shows that the effect of convolution with this PSF is that thin structures become thicker with less contrast. 3.3.2 The edge spread function The ESF (edge spread function) is even easier to measure than the LSF. If half of the image plane is blocked by an Figure 23. The edge spread function of disk function. object, and the other half is not, the resulting image shows a) Input image, b) Output image. the edge spread function. In an X-ray system this is easily The conclusion is that the ESF is the integral of the LSF. done with an aluminium plate. The reverse of integration is differentiating. Thus, a corollary is: The image model of such a situation is a step function. If aligned along the y-axis and located at the origin, the model d is: lsf ( x) = esf ( x) (38) dx 1 if x ≥ 0 f ( x, y ) =  (35) Example: the edge spread function of a disk 0 if x < 0 The LSF for a disk function is given in eq (34). By integration we find the ESF. If x < − R , the LSF is zero, and See Figure 23. Note that this image does not depend on y. thus the ESF is also zero there. Beyond x = − R , the LSF is Therefore, we can also write f ( x) instead of f ( x, y ). positive. Therefore, the ESF will increase until it reaches Consequently, the ESF will not depend y. The response on x = + R. After that, the ESF will be constant. a system with PSF h( x, y ) follows from substitution in eq (26): In the interval − R < x ≤ R , the following integral needs to ∞ ∞ be solved: = esf ( x) ∫ ∫ f ( x − ξ )h(ξ ,η )dη d ξ x x ∫= ∫ ξ = −∞ η = −∞ lsf (ξ )d ξ 2 A R 2 − ξ 2 dξ (39) ∞ ∞ ξ= −R = −R ξ = ∫ ξ = −∞ f (x − ξ ) ∫ η = −∞ h(ξ ,η )dη d ξ (36) ∞ Using a table with integrals or, better, using a symbolic = ∫ ξ = −∞ f ( x − ξ ) lsf (ξ )d ξ mathematical manipulator, e.g. MATLAB’s Symbolic Math Toolbox, the integral can be found. The solution is: The function f ( x − ξ ) is zero if x − ξ < 0 , which is the case x x if ξ ≥ x. Otherwise, the function is 1. Thus: ξ ∫ = −R lsf (ξ )d ξ = 12 π R 2 A + Ax R 2 − x 2 + AR 2 asin   R (40) x esf ( x, y ) = ∫ lsf (ξ )d ξ (37) A further simplification occurs when the volume of the disk ξ = −∞ = function is required to be one: V π= R 2 A 1. In that case, the maximum of the ESF is esf ( R) = 1. Putting it all together:  11   2 2 2 0 if x < − R h =  0 0 0 ,   −2 −2 −2    an input image f , and output image g. Then: esf ( x) 1 if x > + R (41)  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0  1 + x R − x + 1 asin  x  elsewhere = f = 2 2 0 0 0 0 0  gives g  0 0 4 4 4. 2 π R 2 π  R  0 0 0 2 0  0 0 0 0 0  0 0 0 0 0  0 0 −4 −4 −4  The 3 × 3 neighbourhood of the single point with grey level Figure 23 is a graphical representation of the found ESF. 2 has been replaced by the PSF multiplied by 2. If the input The figure shows that the effect of convolution with this PSF consists of two points, like: is that the step becomes smoother, but the step height 0 0 0 0 0  2 2 2 0 0 remains the same. 0 1 0 0 0  0 0 0 0 0 = f = 0 0 0 0 0  gives g  −2 −2 2 4 4  , Homework 2_3: 0 0 0 2 0  0 0 0 0 0 0 0 0 0 0  0 0 −4 −4 −4  Give expressions for the LSF and ESF associated with the rectangular PSF given in eq (9). Assume that the the 3 × 3 neighbourhoods of both are replaced, but since volume is one. the two neighbourhoods overlap, the individual responses Give expressions for the LSF and ESF associated with of the two points are added. the Gaussian PSF given in eq (10). Assume that the As in the continuous case, the discrete convolution is volume is one. = commutative, g n , m f= n , m * hn , m hn , m * f n , m , so that the operation can also be expressed as: 3.4 Discrete convolution +I +J For digital image processing and image analysis, a discrete g n,m = ∑∑h i, j f n −i , m − j (43) version of the convolution integral is needed. The inputs of i= −I j= −J the discrete convolution are: The difference between the two versions of convolution, eq A N × M digital image, represented by a 2D array f n , m (42) and (43), is demonstrated in Figure 24. with n = 1, , N and m = 1, , M. For brevity’s sake, In Figure 24a, the position of the input image is fixed. the digital image is sometimes also denoted by a single The PSF is mirrored and moved to a position (n, m) , symbol f. In MATLAB or Python, notations like f will be and then the calculation as in eq (42) is performed. used for the whole image. To address a specific pixel, In Figure 24b, the position of the PSF is fixed. The input we use: image is mirrored and shifted over a displacement MATLAB: f(m,n) with f n , m  f(m-1,n-1) (n, m) after which the calculation as in eq (43) is Python: f[m,n] with f n , m  f[m,n] performed. A PSF hi , j. For simplicity of notation, we assume that the size of this 2D array is odd: (2 I + 1) × (2 J + 1). The The two versions are mathematical equivalent. elements are enumerated by i = − I , , I and j = − J , , J. Sometimes the notation for the full PSF 3.4.1 Implementation details is also h. MATLAB and Python sample code for convolution The discrete version of the convolution integral in eq (26) Below, sample code is given for image filtering based on becomes: convolution in MATLAB and Python’s OpenCV. n+ I m+ J Listing M 1: Convolution g n,m = ∑ ∑ f p , q hn − p , m − p (42) p =− n I q= m− J %% sample code - convolution close all clear variables This expression is to be evaluated for n = 1, , N and m = 1, , M , so that the output image has size N × M. %% read image and show Iin = imread("..\images\bos.jpg"); To illustrate the operation, suppose that we have a PSF h Iin = im2double(Iin); defined by: figure(1); imshow(Iin)  12 np.flip(hpsf), borderType=cv2.BORDER_REPLICATE) Iout = cv2.normalize(Iout,None,0,1, cv2.NORM_MINMAX) cv2.imshow("Sobel",Iout) cv2.waitKey() Iout = np.uint8(255*Iout) cv2.imwrite('../images/trees_sobel.jpg',Iout) In MATLAB, the core of the code is the function imfilter(). In Python-OpenCV in cv2.filter2D(). Some details are discussed below. - The definition of the point spread function In the example, the PSF hpsf is a 3x3 array, it is convenient to allow negative indices in the 2D PSF array hi , j , where i and j run from − I to + I and from − J to + J. The size is (2 I + 1) × (2 J + 1). By doing so, the origin of the PSF is defined in a natural way. The element h0,0 is automatically situated in the centre of the PSF. However, when implementing the convolution in MATLAB or Python, the indices are not allowed to be negative. Given a MATLAB or Python array hpsf with size JxI, one has to define where the origin is located. If J and I are odd, then by default the origin is in the middle: in MATLAB at (J+1)/2. and (I+1)/2. If they are even, the assumed origin is at L/2 and K/2. In Python, one can change this location by means of a variable called anchor. Figure 24. Two versions of the discrete convolution. - Correlation versus convolution a) The image is fixed; the mirrored PSF moves. By default, the functions imfilter() and cv2.filter2D() b) The PSF is fixed; the mirrored image moves. perform correlation rather than convolution. Correlation is %% create PSF the operator ΣΣ f p.q hn + p , m + q , whereas convolution is the hpsf = [1 0 -1; 2 0 -2; 1 0 -1]; % Sobel kernel operator ΣΣ f p.q hn − p , m − q. Correlation becomes convolution if %% convolve, normalize, show, and write to file one mirrors the PSF across the origin. That is hiconv ,j = h−corr i,− j. Iout = imfilter(Iin,hpsf,'conv','replicate'); In MATLAB, this is done by adding the option ‘conv’ in the Iout = mat2gray(Iout); function imfilter(). In Python’s cv2.filter2D(), this figure(2); imshow(Iout) option does not exist, and one has to mirror the PSF imwrite(Iout,"..\images\trees_sobel.jpg") explicitly with the function np.flip(). Listing P 1: Convolution - Missing data near the image boundary To calculate output pixels that are near the boundary of the # sample code - convolution import cv2 image, input pixels outside the image area are needed. This import numpy as np is illustrated in Figure 25. The pixels outside the image area are not available. To circumvent the problem, these outside #%% read and show image Iin = cv2.imread('../images/trees.jpg') pixels are extrapolated. For that, different options exist. The Iin = Iin/255. most important ones are: cv2.imshow("trees",Iin) Reflection: Pixels outside the bounds of the image #%% create Sobel PSF area are computed by mirror-reflecting the array across hpsf = np.array([[1,0,-1],[2,0,-2],[1,0,-1]]) the image border. MATLAB: ‘symmetric’, Python OpenCV: cv2.BORDER_REFLECT. #%% convolve, normalize, show, and write to file Iout = cv2.filter2D(Iin,-1,  13 The proof is easily given by substitution of the expression for convolution. Cascading convolutions Suppose that an image f is convolved with a PSF h1 after which the result is convolved once again with a PSF h 2 : g = ( f * h1 ) * h 2 (45) One can perform this cascaded convolution in one step by first calculating a new PSF h = h1 * h 2. In other words: Figure 25. Missing data near the boundary of the image. =g (= f * h1 ) * h 2 f *(h1 * h 2 ) (46) Replication: Pixels outside the bounds of the image The proof is easily given in the frequency domain using the area are assumed to be equal to the nearest pixel property stated in eq (85). See Section 4.2. inside the image area. MATLAB: ‘replicate’, Python’s OpenCV: cv2.BORDER_REPLICATE. Separability Constant: Pixels outside the bounds of the image area If the PSF is separable, e.g. hi , j = ai b j , then the 2D are given a defined value. In MATLAB: just give the convolution can be simplified to two 1D convolutions: value as third argument in the function call, Python’s +I +J OpenCV: add the option cv2.BORDER_CONSTANT. g n,m = ∑ ∑ ab i= −I j= −J i j f n −i , m − j +I +J (47) Homework 2_4: = Create a MATLAB file H2_4.m that reads and shows the file ∑ = i= aw −I with w ∑b i n −i , m n,m j= −J j n,m − j f checker_board.png. Create a PSF with the following code: This reduces the computational complexity enormously. hpsf = ones(11)/11^2; Computational complexity The PSF is a discrete rectangular function. Effectively, the Suppose that we have a PSF of size (2 I + 1) × (2 J + 1) which convolution with this PSF performs averaging of the pixels in the MATLAB or Python realization becomes a JxI PSF. in each 11x11 neighbourhood. Apply this filter with: with = I 2 I + 1 and= J 2 J + 1. The convolution the option ‘replicate’, the option ‘symmetric’, +I +J and with no option of such kind. g n,m = ∑∑h i= −I J = −J i, j f n −i , m − j Save the resulting images to file in folder..\images. Observe and explain the differences. Hint: read the MATLAB consists of an inner loop in which J operations have to be documentation to see what the default behaviour is when accomplished. An operation consists of a multiplication and no option is given. Note: don’t forget to add the option an addition. This is repeated I times in the second loop. In ‘conv’ in all cases as we want convolution and not full, there are IJ operations. However, this is only for one correlation. output pixel g n , m. As the output image has size MxN, in total, the number of operations is NMIJ. 3.4.2 Properties of convolution As a numerical example, suppose that a machine can Linearity of convolution perform 1 operation in 0.1 ns , then a convolution on an Suppose that an image f is formed as a linear image of size 1500x1500 with a PSF of 11x11 will take combination of two other images f1 and f 2. Thus, around 0.03 s. f α f1 + β f 2. The image is convolved with a PSF h. = According to the linearity property of convolution: If the PSF is separable, eq (47) applies. The first convolution takes NMJ operations; the second one NMI. In (α f1 + β f 2 ) * h = α (f1* h) + β (f 2* h) (44) full this is NM(I+J). In the example, the time needed for the convolution reduces from 0.03 s to 0.005 s.  14 a) b) c) a) b) Figure 26. Noise reduction = low-pass filtering. Figure 27. High pass filtering. a) Original, b) smoothed; R=2, c) smoothed; R=3. a) R=2, b) R=3. a) b) 3.5 Applications The two most important applications of convolution are image enhancement and feature extraction for image analysis. 3.5.1 Image enhancement and noise reduction Noise reduction Figure 28. Image sharpening with α = 2. Real images are always subject to noise. In X-ray images, a) R=2, b) R=3. for instance, the image is built up from a finite number of X- textures, especially in regions without contrast. Application ray quanta which causes quantum noise. A PSF is essential of a convolution with a disk-like PSF, with two different in the formation of the image to reduce the noise. With radii, are shown in Figure 26b and c. The quantum noise discrete convolution, the spatial extent of the PSF can be decreases as the radius increases. This is at the cost of controlled. more blurry lines and edges due to the LSF and ESF associated with the PSF. Convolutions with the PSFs mentioned in Section 2.3 all have the effect to reduce the noise. The requirement is that For reasons that become clear in the next sections filtering the mean grey value over the whole image will not be with the mentioned PSFs is called low-pass filtering. changed by the convolution. This is achieved by choosing the amplitude of the PSF such that its volume becomes High-pass filtering and image sharpening one. For instance, in case of a disk function the amplitude In physical imaging systems, it is difficult to create PSFs must be A = 1 (π R 2 ). The convolution then effectively with negative values, but in digital domain it is easy. This averages the grey level in a circular neighbourhood with opens possibilities for many applications. One is image radius R. Similar behaviour is seen in the rectangular PSF. sharpening. The Gaussian PSF, with A = 1 (2πσ 2 ) , also implements To introduce this, first consider high-pass filtering, which is averaging, However, this is a weighed averaging. Pixels the complement of low-pass filtering: further away from the central pixels are weighed less. In digital image processing, the PSFs must be discretized. g hpf = f − g lpf For a disk-like PSF with radius R = 2 , the PSF becomes a = f − f * hlpf (49) 5x5 array: = f * 1 − f * hlpf  0 0.2 0.5 0.2 0 = f * h hpf with h hpf = 1 − hlpf  0.2 1 1 1 0.2   0.5 1 0.5  1 h disk = 1 1 (48) The PSF 1 is the identity PSF which leaves the image 12.6  0.2 1 0.2   0 1 1 unchanged. That is, for instance:  0.2 0.5 0.2 0  Pixels that are near the edge of the disk are interpolated. It 0 0 0 0 0 0 0 0 0 0 can be seen that ∑∑ h disk = 1. 1 = 0 0 0 1 0 0. (50) 0 0 0 0 0 0 Figure 26a is a snippet from a fluoroscopy image. Upon  0 0 0  inspection, the effect of quantum noise is seen as woolly  15 hlpf is a PSF that realizes low-pass filtering. Examples are the discrete versions of the PSFs in Section 2.3, which can also be used for noise reduction. By application of eq (49) and (50) in the previous example, the disk-like PSF with R = 2 turns into a high-pass filter as follows:  0 −0.2 −0.5 −0.2 0  −0.2 −1 −1 −1 −0.2   −0.5 −1 −0.5  1 h hpf = −1 11.6 (51) 12.6  −0.2 −1 −1 −1 −0.2   0 0  −0.2 −0.5 −0.2  Figure 29. Sobel operator: vertical edges. The resulting images, are shown in Figure 27, High-pass Homework 2_5: creation of PSFs filtering extracts the fine detailed structures in an image, In MATLAB, create a file H2_5.m. Initialize it properly. Then: and suppresses all slowly fluctuating structures. create a disk-like PSF hdisk with radius R=10 with the function fspecial(). Visualize it as a 3D surface mesh. To use this for image sharpening, the following operation For that, use the following code: can be used: %% create psfs f + α g hpf = g sharp = f + α f * h hpf (52) hdisk = fspecial('disk',10); The operation amplifies the detailed structures in the %% create a figure window with defined size image without changing the slowly fluctuation structures. f = figure(1); % get figure handle f.Position = [200,200,1000,300]; % window size The parameter α is the amplification factor. Examples are shown in Figure 28. %% divide the window area in 3 tiles %% and select the first one 3.5.2 Image feature extraction subplot(1,3,1) Another application of convolution is the extraction of %% show psf as surface mesh mesh(hdisk,"EdgeColor","black") specific image features. One of the older examples is the xlabel('x'); ylabel('y') Sobel operator, which aims to extract the edges of objects title('disk') in the image. There is a PSF for vertical edges and one for Extend the code to generate two other PSFs. The function horizontal edges: fspecial() generates predefined PSFs. Use this function:  1 2 1  1 0 −1 to create a Gaussian PSF, called hgauss. The size h horz = 0 0 0 h vert  2 0 −2  (53) should be 41x41, and the sigma must be 5. Read the  −1 −2 −1  1 0 −1  documentation. Visualize this PSF in the second tile of the figure window. Application of the vertical Sobel operator is shown in Figure to create a rectangular window, called hrect. The size 29. must be 18x18. Use fspecial(). Visualize the PSF in These types of convolutions are fundamental for image the third tile. analysis. Operations for finding edges, lines, corners, etc in Print the figure to file. For that, you can use: print('-r300','-dpng','..\images\H2_5_PSFs.png') images are most rigorously treated in the so-called scale- space theory. This theory is also the basis for key point Apply the 3 PSFs to the test image impoints.png. Don’t detection: detection of specific points in an image. forget to use the function im2double(). (Why?) Show the resulting 3 images to 3 tiles in a figure window in the same Another branch of image analysis is based on deep way as the 3 meshes. Apply grey level windowing to make learning. The symbiosis of neural networks and the sure that the results are visible. Print the figure window to convolution operator has led to so-called convolutional file. neural networks, CNNs. Many image processing tasks, such as segmentation, image classification, image registration, Observe the results and explain. image generation, are using these CNNs.  16 Homework 2_6: noise reduction PSF with fixed σ = 5. We treat α as a design parameter The file test_image.png contains four subimages. Two that has to be chosen yet. subimages are test pattern consisting of a concentric wave The exercise is to be implemented in Python. OpenCV does function and a radial wave function. A third image consists not have a function to generate general-purpose PSFs. of a disk and rings. The fourth image is a bronchoscopic Python code to build a Gaussian 1 PSF from scratch is: image. #%% create a PSF A second file, test_image_noisy.png, is the same as the I=20 clean test image except that it has been contaminated by s = 4. i = np.linspace(-I,I,2*I+1) noise. The goal of this exercise is to assess which of the xy = np.meshgrid(i,i) PSFs created in Homework 2_5 is best in reducing the R2 = xy**2 + xy**2 noise and at the same time is preserving the detailedness hgauss = np.exp(-R2/(2*s**2))/(2*np.pi*s**2) in the image as much as possible. The function meshgrid() creates two 2D arrays with the The first assessment will be based on visible inspection. coordinates of the grid points of the PSF array. Such a human visual perception is subjective. A second, The design is tweaked with the parameter alpha. To have a more objective assessment is based on the interpretation more objective criterion for selecting alpha, plots of the of a plot of a cross-section of the image data over an cross-sections of the blurred image and the restored image interesting part of the test image. As an example, to create can be compared. such a plot for the clean image, you can use: To plot the cross-sections, simply use plt.plot() from the im = imread('..\images\test_image.png'); matplotlib package. For that you need to import the im = im2double(im); xsection = im(376,376:1068);% get cross-section subpackage: figure; import matplotlib.pyplot as plt plot(xsection); % plot data To save a plot to file, and to activate displaying the plot in a Read the image test_image_noisy.png. Apply figure window, use (in that order): convolution with the PSF hdisk created in Homework plt.savefig("../images/H2_7plt.png") 2_5. Write the image to file in..\images. Extract the plt.show() cross-section as was done in the clean image. Repeat the procedure, but now with the PSF hgauss, Implement the image sharpening and the plotting of and the PSF hrect. the cross-section. Apply to testimage_blurry.png. Plot the four cross-sections in one figure. Print the Initially, use σ = 5 and α = 1. figure to file in the folder..\images. Note, to plot Variate the parameter α. By inspecting the resulting multiple graphs in one axis, call plot of the cross-section, choose the value that suits hold on you best. Which criterion(s) do you use? after the first call to plot(). Write the final image to../images/H2_7.jpg. Don’t forget to add labels, a title, or a legend with Observe and explain all relevant aspects that you see. xlabel(), ylabel(), title() or legend(). By visually inspecting the plot with the four graphs, rank the three PSFs according to their ability to reduce noise. If the PSFs perform equally well, then say so. By visually inspecting the plot with the four graphs, rank the PSFs according to their ability to preserve sharpness. If PSFs perform equally well, then say so. Observe, mention and try to explain all remarkable things that you see in the results. Homework 2_7: image sharpening The file test_image_blurry.png is a blurred version of the 1 By the way, Python code to generate a disk-like PSF is: original image in test_image.png. The goal of this exercise hdisk = np.float64(R2