Question? Leave a message!




Introduction to Image Processing

Introduction to Image Processing 8
Introduction to Image Processing Prof. George Wolberg Dept. of Computer Science City College of New YorkCourse Description •Intense introduction to image processing. •Intended for advanced undergraduate and graduate students. •Topics include: Image enhancement Digital filtering theory, Fourier transforms Image reconstruction, resampling, antialiasing Scanline algorithms, geometric transforms Warping, morphing, and visual effects 2 Wolberg: Image Processing Course NotesSyllabus Topic Week Introduction / overview 1 Point operations 23 Neighborhood operations 4 Fast Fourier transforms (FFT) 56 Sampling theory 78 Midterm, Image reconstruction 9 Fast filtering for resampling 10 Spatial transformations, texture mapping 11 Separable warping algorithms; visual effects 1214 3 Wolberg: Image Processing Course NotesRequired Text Rafael Gonzalez and Richard Woods, Digital rd Image Processing, 3 Edition, Prentice Hall, Wesley, 2008. 4 Wolberg: Image Processing Course NotesSupplementary Texts Milan Sonka, Vaclav Hlavac, and Roger Boyle, Image Processing, Analysis, and Machine Vision, Cengage Learning, 2014. George Wolberg, Digital Image Warping, IEEE Computer Society Press, 1990. 5 Wolberg: Image Processing Course NotesGrading •The final grade is computed as follows: Midterm exam: 25 Final exam: 25 Homework programming assignments: 50 •Substantial programming assignments are due every three weeks. •Proficiency in C/C++ is expected. •Prereqs: CSc 22100 6 Wolberg: Image Processing Course NotesContact Information •Prof. Wolberg Office hours: After class and by appointment Email: wolbergcs.ccny.cuny.edu •Teaching Assistant (TA): Siavash Zokai Email: ccny.cs470gmail.com •See class web page for all class info such as homework and sample source code: wwwcs.ccny.cuny.edu/wolberg/cs470 7 Wolberg: Image Processing Course NotesObjectives •These notes accompany the textbooks: “Digital Image Processing” by Gonzalez/Woods “Digital Image Warping” by George Wolberg •They form the basis for approximately 14 weeks of lectures. •Programs in C/C++ will be assigned to reinforce understanding of the material. Four homework assignments Each due in 3 weeks and requiring 4 programs 8 Wolberg: Image Processing Course NotesWhat is Image Processing Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives •In this lecture we: Explore what image processing is about Compare it against related fields Provide historical introduction Survey some application areas 10 Wolberg: Image Processing Course NotesWhat is Digital Image Processing • Computer manipulation of pictures, or images, that have been converted into numeric form. Typical operations include: Contrast enhancement Remove blur from an image Smooth out graininess, speckle, or noise Magnify, minify, or rotate an image (image warping) Geometric correction Image compression for efficient storage/transmission 11 Wolberg: Image Processing Course NotesImage Processing Goals • Image processing is a subclass of signal processing concerned specifically with pictures • It aims to improve image quality for human perception: subjective computer interpretation: objective • Compress images for efficient storage/transmission 12 Wolberg: Image Processing Course NotesRelated Fields Image Processing Image Computer Computer Graphics Vision Scene Description 13 Wolberg: Image Processing Course NotesOverlap with Related Fields Image Processing Image Noise reduction Contrast enhancement Filtering Texture mapping Lowlevel Imagein / Imageout Antialiasing Extract attributes Computer Computer Edge detection Midlevel Graphics Vision Segmentation Imagein / Featuresout Recognition Highlevel Cognitive functions Scene Description 14 Wolberg: Image Processing Course NotesDistinctions • No clear cut boundaries between image processing on the one end and computer vision at the other • Defining image processing as imagein/imageout does not account for computation of average intensity: imagein / numberout image compression: imagein / coefficientsout • Nevertheless, imagein / imageout is true most of time Output Image Description Input Image Computer Image Processing Vision Computer Artificial Description Graphics Intelligence 15 Wolberg: Image Processing Course NotesImage Processing: 19601970 Geometric correction and image enhancement applied to Ranger 7 pictures of the moon. Work conducted at the Jet Propulsion Laboratory. 16 Wolberg: Image Processing Course NotesImage Processing: 19701980 • Invention of computerized axial tomography (CAT) • Emergence of medical imaging • Rapid growth of Xray imaging for CAT scans, inspection, and astronomy • LANDSAT earth observation 17 Wolberg: Image Processing Course NotesImage Processing: 19801990 • Satellite infrared imaging: LANDSAT, NOAA • Fast resampling and texture mapping 18 Wolberg: Image Processing Course NotesImage Processing: 19902000 • Morphing / visual effects algorithms • JPEG/MPEG compression, wavelet transforms • Adobe PhotoShop 19 Wolberg: Image Processing Course NotesImage Processing: 2000 • Widespread proliferation of fast graphics processing units (GPU) from nVidia and ATI to perform realtime image processing • Ubiquitous digital cameras, camcorders, and cell phone cameras rely heavily on image processing and compression 20 Wolberg: Image Processing Course NotesSources of Images • The principal energy source for images is the electromagnetic energy spectrum. • EM waves = stream of massless (proton) particles, each traveling in a wavelike pattern at the speed of light. Spectral bands are grouped by energy/photon Gamma rays, Xrays, UV, Visible, Infrared, Microwaves, radio waves • Other sources: acoustic, ultrasonic, electronic 21 Wolberg: Image Processing Course NotesGammaRay Imaging • Used in nuclear medicine, astronomy • Nuclear medicine: patient is injected with radioactive isotope that emits gamma rays as it decays. Images are produced from emissions collected by detectors. 22 Wolberg: Image Processing Course NotesXRay Imaging • Oldest source of EM radiation for imaging • Used for CAT scans • Used for angiograms where Xray contrast medium is injected through catheter to enhance contrast at site to be studied. • Industrial inspection 23 Wolberg: Image Processing Course NotesUltraviolet Imaging • Used for lithography, industrial inspection, flourescence microscopy, lasers, biological imaging, and astronomy • Photon of UV light collides with electron of fluorescent material to elevate its energy. Then, its energy falls and it emits red light. 24 Wolberg: Image Processing Course NotesVisible and Infrared Imaging (1) • Used for astronomy, light microscopy, remote sensing 25 Wolberg: Image Processing Course NotesVisible and Infrared Imaging (2) • Industrial inspection inspect for missing parts missing pills unacceptable bottle fill unacceptable air pockets anomalies in cereal color incorrectly manufactured replacement lens for eyes 26 Wolberg: Image Processing Course NotesMicrowave Imaging • Radar is dominant application • Microwave pulses are sent out to illuminate scene • Antenna receives reflected microwave energy 27 Wolberg: Image Processing Course NotesRadioBand Imaging • Magnetic resonance imaging (MRI): places patient in powerful magnet passes radio waves through body in short pulses each pulse causes a responding pulse of radio waves to be emitted by patient’s tissues Location and strength of signal is recorded to form image 28 Wolberg: Image Processing Course NotesImages Covering EM Spectrum 29 Wolberg: Image Processing Course NotesNonEM modality: Ultrasound • Used in geological exploration, industry, medicine: transmit highfreq (15 MHz) sound pulses into body record reflected waves calculate distance from probe to tissue/organ using the speed of sound (1540 m/s) and time of echo’s return display distance and intensities of echoes as a 2D image 30 Wolberg: Image Processing Course NotesNonEM modality: Scanning Electron Microscope • Stream of electrons is accelerated toward specimen using a positive electrical potential • Stream is focused using metal apertures and magnetic lenses into a thin beam • Scan beam; record interaction of beam and sample at each location (dot on phosphor screen) 31 Wolberg: Image Processing Course NotesVisible Spectrum • Thin slice of the full electromagnetic spectrum 32 Wolberg: Image Processing Course NotesHuman Visual System Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives •In this lecture we discuss: Structure of human eye Mechanics of human visual system (HVS) Brightness adaptation and discrimination Perceived brightness and simultaneous contrast 34 Wolberg: Image Processing Course NotesHuman and Computer Vision • We observe and evaluate images with our visual system • We must therefore understand the functioning of the human visual system and its capabilities for brightness adaptation and discrimination: What intensity differences can we distinguish What is the spatial resolution of our eye How accurately do we estimate distances and areas How do we sense colors By which features can we detect/distinguish objects 35 Wolberg: Image Processing Course NotesExamples Parallel lines : 5 variation in length Circles: 10 variation in radii Upper line falsely appears longer Vertical line falsely appears longer 36 Wolberg: Image Processing Course NotesStructure of the Human Eye • Shape is nearly spherical • Average diameter = 20mm • Three membranes: Cornea and Sclera Choroid Retina 37 Wolberg: Image Processing Course NotesStructure of the Human Eye: Cornea and Sclera • Cornea Tough, transparent tissue that covers the anterior surface of the eye • Sclera Opaque membrane that encloses the remainder of the optical globe 38 Wolberg: Image Processing Course NotesStructure of the Human Eye: Choroid • Choroid Lies below the sclera Contains network of blood vessels that serve as the major source of nutrition to the eye. Choroid coat is heavily pigmented and hence helps to reduce the amount of extraneous light entering the eye and the backscatter within the optical globe 39 Wolberg: Image Processing Course NotesLens and Retina • Lens Both infrared and ultraviolet light are absorbed appreciably by proteins within the lens structure and, in excessive amounts, can cause damage to the eye • Retina Innermost membrane of the eye which lines the inside of the wall’s entire posterior portion. When the eye is properly focused, light from an object outside the eye Is imaged on the retina. 40 Wolberg: Image Processing Course NotesReceptors • Two classes of light receptors on retina: cones and rods • Cones 67 million cones lie in central portion of the retina, called the fovea. Highly sensitive to color and bright light. Resolve fine detail since each is connected to its own nerve end. Cone vision is called photopic or brightlight vision. • Rods 75150 million rods distributed over the retina surface. Reduced amount of detail discernable since several rods are connected to a single nerve end. Serves to give a general, overall picture of the field of view. Sensitive to low levels of illumination. Rod vision is called scotopic or dimlight vision. 41 Wolberg: Image Processing Course NotesDistribution of Cones and Rods • Blind spot: no receptors in region of emergence of optic nerve. • Distribution of receptors is radially symmetric about the fovea. • Cones are most dense in the center of the retina (e.g., fovea) • Rods increase in density from the center out to 20° and then decrease 42 Wolberg: Image Processing Course NotesBrightness Adaptation (1) • The eye’s ability to discriminate between intensities is important. • Experimental evidence suggests that subjective brightness (perceived) is a logarithmic function of light incident on eye. Notice approximately linear response in logscale below. Wide range of intensity levels to which HVS can adapt: from scotopic threshold to glare limit (on the order of 1010) Range of subjective brightness that eye can perceive when adapted to level B a 43 Wolberg: Image Processing Course NotesBrightness Adaptation (2) • Essential point: the HVS cannot operate over such a large range simultaneously. • It accomplishes this large variation by changes in its overall sensitivity: brightness adaptation. • The total range of distinct intensity levels it can discriminate simultaneously is rather small when compared with the total adaptation range. • For any given set of conditions, the current sensitivity level of the HVS is called the brightness adaptation level (B in figure). a 44 Wolberg: Image Processing Course NotesBrightness Discrimination (1) • The ability of the eye to discriminate between intensity changes at any adaptation level is of considerable interest. • Let I be the intensity of a large uniform area that covers the entire field of view. • Let ΔI be the change in object brightness required to just distinguish the object from the background. • Good brightness discrimination: ΔI / I is small. • Bad brightness discrimination: ΔI / I is large. • ΔI / I is called Weber’s ratio. 45 Wolberg: Image Processing Course NotesBrightness Discrimination (2) • Brightness discrimination is poor at low levels of illumination, where vision is carried out by rods. Notice Weber’s ratio is large. • Brightness discrimination improves at high levels of illumination, where vision is carried out by cones. Notice Weber’s ratio is small. rods cones 46 Wolberg: Image Processing Course NotesChoice of Grayscales (1) • Let I take on 256 different intensities: 0  I 1 for j = 0,1,…,255. j • Which levels we use Use eye characteristics: sensitive to ratios of intensity levels rather than to absolute values (Weber’s law: B/B = constant) For example, we perceive intensities .10 and .11 as differing just as much as intensities .50 and .55. 47 Wolberg: Image Processing Course NotesChoice of Grayscales (2) • Levels should be spaced logarithmically rather than linearly to achieve equal steps in brightness: 2 3 255 I , I = rI , I = rI = r I , I = rI = r I , …., I = r I =1, 0 1 0 2 1 0 3 2 0 255 0 where I is the lowest attainable intensity. 0 1/255 • r = (1/I ) , 0 j j/255 (1j/255) (255j)/255 • I = r I = (1/I ) I = I = I j 0 0 0 0 0 • In general, for n+1 intensities: 1/n, r = (1/I ) 0 (nj)/n I = I for 0 jn j 0 48 Wolberg: Image Processing Course NotesChoice of Grayscales (3) • Example: let n=3 and I =1/8: 0 r = 2 (3/3) I = (1/8) 0 (2/3) I = (1/8) = ¼ 1 (1/3) I = (1/8) = ½ 2 (0/3) I = (1/8) =1 3 For CRTs, 1/200 I 1/40. 0 I0 because of light reflection from the phosphor 0 within the CRT. Linear grayscale is close to logarithmic for large number of graylevels (256). 49 Wolberg: Image Processing Course NotesPerceived Brightness • Perceived brightness is not a simple function of intensity. • The HVS tends to over/undershoot around intensity discontinuities. • The scalloped brightness bands shown below are called Mach bands, after Ernst Mach who described this phenomenon in 1865. 50 Wolberg: Image Processing Course NotesSimultaneous Contrast (1) • A region’s perceived brightness does not depend simply on its intensity. It is also related to the surrounding background. 51 Wolberg: Image Processing Course NotesSimultaneous Contrast (2) • An example with colored squares. 52 Wolberg: Image Processing Course NotesProjectors •Why are projection screens white Reflects all colors equally well •Since projected light cannot be negative, how are black areas produced Exploit simultaneous contrast The bright area surrounding a dimly lit point makes that point appear darker 53 Wolberg: Image Processing Course NotesVisual Illusions (1) 54 Wolberg: Image Processing Course NotesVisual Illusions (2) • Rotating snake illusion • Rotation occurs in relation to eye movement • Effect vanishes on steady fixation • Illusion does not depend on color • Rotation direction depends on the polarity of the luminance steps • Asymmetric luminance steps are required to trigger motion detectors 55 Wolberg: Image Processing Course NotesDigital Image Fundamentals Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives •In this lecture we discuss: Image acquisition Sampling and quantization Spatial and graylevel resoIution 57 Wolberg: Image Processing Course NotesSensor Arrangements • Three principal sensor arrangements: Single, line, and array 58 Wolberg: Image Processing Course NotesSingle Sensor • Photodiode: constructed of silicon materials whose output voltage waveform is proportional to light. • To generate a 2D image using a single sensor, there must be relative displacements in the horizontal and vertical directions between the sensor and the area to be imaged. • Microdensitometers: mechanical digitizers that use a flat bed with the single sensor moving in two linear directions. 59 Wolberg: Image Processing Course NotesSensor Strips • Inline arrangement of sensors in the form of a sensor strip. • The strip provides imaging elements in one direction. • Motion perpendicular to strip images in the other direction. • Used in flat bed scanners, with 4000 or more inline sensors. 60 Wolberg: Image Processing Course NotesSensor Arrays • Individual sensors are arranged in the form of a 2D array. • Used in digital cameras and camcorders. • Entire image formed at once; no motion necessary. 61 Wolberg: Image Processing Course NotesSignals • A signal is a function that conveys information 1D signal: f(x) waveform 2D signal: f(x,y) image 3D signal: f(x,y,z) volumetric data or f(x,y,t) animation (spatiotemporal volume) 4D signal: f(x,y,z,t) snapshots of volumetric data over time • The dimension of the signal is equal to its number of indices. • In this course, we focus on 2D images: f(x,y) • Efficient implementation often calls for 1D row or column processing. That is, process the rows independently and then process the columns of the resulting image. 62 Wolberg: Image Processing Course NotesDigital Image • Image produced as an array (the raster) of picture elements (pixels or pels) in the frame buffer. 63 Wolberg: Image Processing Course NotesImage Classification (1) • Images can be classified by whether they are defined over all points in the spatial domain and whether their image values have finite or infinite precision. • If the position variables (x,y) are continuous, then the function is defined over all points in the spatial domain. • If (x,y) is discrete, then the function can be sampled at only a finite set of points, i.e., the set of integers. • The value that the function returns can also be classified by its precision, independently of x and y. 64 Wolberg: Image Processing Course NotesImage Classification (2) • Quantization refers to the mapping of real numbers onto a finite set: a manytoone mapping. • Akin to casting from double precision to an integer. Space Image Values Classification continuous continuous analog (continuous) image continuous discrete intensity quantization discrete continuous spatial quantization discrete discrete digital (discrete) image 65 Wolberg: Image Processing Course NotesImage Formation • The values of an acquired image are always positive. There are no negative intensities: 0 f(x,y) ∞ • Continuous function f(x,y) = i(x,y) r(x,y), where 0 i(x,y) ∞ is the illumination 0 r(x,y) 1 is the reflectance of the object • r(x,y) = 0 is total absorption and r(x,y) = 1 is total reflectance. • Replace r(x,y) with transmissivity term t(x,y) for chest Xray. 66 Wolberg: Image Processing Course NotesTypical Values of Illumination and Reflectance 2 • The following i(x,y) illumination values are typical (in lumens/m ): Illumination of sun on Earth on a clear day: 90,000 Illumination of sun on Earth on a cloudy day: 10,000 Illumination of moon on Earth on a clear night: 0.1 Illumination level in a commercial office: 1000 Illumination level of video projectors: 10001500 • The following r(x,y) reflectance values are typical: Black velvet: 0.01 Stainless steel: 0.65 Flatwhite wall paint: 0.80 Silverplated metal: 0.90 Snow: 0.93 67 Wolberg: Image Processing Course NotesGraylevels • The intensity of a monochrome image at any coordinate (x,y) is called graylevel L, where Lmin ≤ L ≤ Lmax • The office illumination example indicates that we may expect Lmin ≈ .01 1000 = 10 (virtual black) Lmax ≈ 1 1000 = 1000 (white) • Interval Lmin, Lmax is called the grayscale. • In practice, the interval is shifted to the 0, 255 range so that intensity can be represented in one byte (unsigned char). • 0 is black, 255 is white, and all intermediate values are different shades of gray varying from black to white. 68 Wolberg: Image Processing Course NotesGenerating a Digital Image • Sample and quantize continuous input image. 69 Wolberg: Image Processing Course NotesImage Sampling and Quantization • Sampling: digitize (discretize) spatial coordinate (x,y) • Quantization: digitize intensity level L 70 Wolberg: Image Processing Course NotesEffects of Varying Sampling Rate (1) • Subsampling was performed by dropping rows and columns. • The number of gray levels was kept constant at 256. 71 Wolberg: Image Processing Course NotesEffects of Varying Sampling Rate (2) • Size differences make it difficult to see effects of subsampling. 72 Wolberg: Image Processing Course NotesSpatial Resolution • Defined as the smallest discernable detail in an image. • Widely used definition: smallest number of discernable line pairs per unit distance (100 line pairs/millimeter). • A line pair consists of one line and its adjacent space. • When an actual measure of physical resolution is not necessary, it is common to refer to an MxN image as having spatial resolution of MxN pixels. 73 Wolberg: Image Processing Course NotesGraylevel Resolution • Defined as the smallest discernable change in graylevel. • Highly subjective process. • The number of graylevels is usually a power of two: k k bits of resolution yields 2 graylevels. When k=8, there are 256 graylevels ← most typical case • Blackandwhite television uses k=6, or 64 graylevels. 74 Wolberg: Image Processing Course NotesEffects of Varying Graylevels (1) • Number of graylevels reduced by dropping bits from k=8 to k=1 • Spatial resolution remains constant. 75 Wolberg: Image Processing Course NotesEffects of Varying Graylevels (2) • Notice false contouring in coarsely quantized images. • Appear as fine ridgelike structures in areas of smooth gray levels. 76 Wolberg: Image Processing Course NotesStorage Requirements • Consider an NxN image having k bits per pixel. • Color (RGB) images require three times the storage (assuming no compression). 77 Wolberg: Image Processing Course NotesLarge Space of Images • Any image can be downsampled and represented in a few bits/pixel for use on small coarse displays. • How many unique images can be displayed on an NxN kbit display k 2 possible values at each pixel 2 N pixels 2 k N Total: (2 ) • This total is huge even for k=1 and N=8: 64 18,446,744,073,709,551,616 ← 2 • It would take 19,498,080,578 years to view this if it were laid out on video at 30 frames/sec. 78 Wolberg: Image Processing Course NotesPoint Operations Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we describe point operations commonly used in image processing: Thresholding Quantization (aka posterization) Gamma correction Contrast/brightness manipulation Histogram equalization/matching 80 Wolberg: Image Processing Course NotesPoint Operations • Output pixels are a function of only one input point: g(x,y) = Tf(x,y) • Transformation T is implemented with a lookup table: An input value indexes into a table and the data stored there is copied to the corresponding output position. The LUT for an 8bit image has 256 entries. g(x,y)=Tf(x,y) LUT Input: f(x,y) Output: g(x,y) 81 Wolberg: Image Processing Course NotesGraylevel Transformation T Contrast enhancement: Thresholding: Darkens levels below m Replace values below m to black (0) Brightens levels above m Replace values above m to white (255) 82 Wolberg: Image Processing Course NotesLookup Table: Threshold • Init LUT with samples taken from thresholding function T LUT 0 0 1 0 2 0 … 0 … 0 m 255 255 … g(x,y)=Tf(x,y) 255 … 255 255 255 Input: f(x,y) Output: g(x,y) 83 Wolberg: Image Processing Course NotesLookup Table: Quantization • Init LUT with samples taken from quantization function T LUT 0 0 … 0 64 64 … 64 128 128 … 128 192 192 g(x,y)=Tf(x,y) 192 255 255 Input: f(x,y) Output: g(x,y) 84 Wolberg: Image Processing Course NotesThreshold Program • Straightforward implementation: // iterate over all pixels for(i=0; itotal; i++) if(ini thr) outi = BLACK; else outi = WHITE; • Better approach: exploit LUT to avoid total comparisons: // init lookup tables for(i=0; ithr; i++) luti = BLACK; for(; iMXGRAY; i++) luti = WHITE; // iterate over all pixels for(i=0; itotal; i++) outi = lutini; 85 Wolberg: Image Processing Course NotesQuantization Program • Straightforward implementation: // iterate over all pixels scale = MXGRAY / levels; for(i=0; itotal; i++) outi = scale (int) (ini/scale); • Better approach: exploit LUT to avoid total mults/divisions: // init lookup tables scale = MXGRAY / levels; for(i=0; iMXGRAY; i++) luti = scale (int) (i/scale); // iterate over all pixels for(i=0; itotal; i++) outi = lutini; 86 Wolberg: Image Processing Course NotesQuantization Artifacts • The false contours associated with quantization are most noticeable in smooth areas. • They are obscured in highly textured regions. Original image Quantized to 8 levels 87 Wolberg: Image Processing Course NotesDither Signal • Reduce quantization error by adding uniformly distributed white noise (dither signal) to the input image prior to quantization. • Dither hides objectional artifacts. • To each pixel of the image, add a random number in the range m, m, where m is MXGRAY/quantizationlevels. v out Uniform noise 0 255 v thr in 3 bpp (8 levels) 8 bpp (256 levels) 88 Wolberg: Image Processing Course NotesComparison 1 bpp 2 bpp 3 bpp 4 bpp Quantization Dither/ Quantization 89 Wolberg: Image Processing Course NotesEnhancement • Point operations are used to enhance an image. • Processed image should be more suitable than the original image for a specific application. • Suitability is applicationdependent. • A method which is quite useful for enhancing one image may not necessarily be the best approach for enhancing another image. • Very subjective 90 Wolberg: Image Processing Course NotesTwo Enhancement Domains • Spatial Domain: (image plane) Techniques are based on direct manipulation of pixels in an image • Frequency Domain: Techniques are based on modifying the Fourier transform of an image • There are some enhancement techniques based on various combinations of methods from these two categories. 91 Wolberg: Image Processing Course NotesEnhanced Images • For human vision The visual evaluation of image quality is a highly subjective process. It is hard to standardize the definition of a good image. • For machine perception The evaluation task is easier. A good image is one which gives the best machine recognition results. • A certain amount of trial and error usually is required before a particular image enhancement approach is selected. 92 Wolberg: Image Processing Course NotesThree Basic Graylevel Transformation Functions • Linear function Negative and identity transformations • Logarithmic function Log and inverselog transformations • Powerlaw function th th n power and n root transformations 93 Wolberg: Image Processing Course NotesImage Negatives • Negative transformation : s = (L–1) – r • Reverses the intensity levels of an image. • Suitable for enhancing white or gray detail embedded in dark regions of an image, especially when black area is large. 94 Wolberg: Image Processing Course NotesLog Transformations s = c log (1+r) • c is constant and r  0 • Log curve maps a narrow range of low graylevels in input image into a wider range of output levels. • Expands range of dark image pixels while shrinking bright range. • Inverse log expands range of bright image pixels while shrinking dark range. 95 Wolberg: Image Processing Course NotesExample of Logarithm Image 6 • Fourier spectrum image can have intensity range from 0 to 10 or higher. • Log transform lets us see the detail dominated by large intensity peak. 6 Must now display 0,6 range instead of 0,10 range. Rescale 0,6 to the 0,255 range. 96 Wolberg: Image Processing Course NotesPowerLaw Transformations  s = cr • c and  are positive constants • Powerlaw curves with fractional values of  map a narrow range of dark input values into a wider range of output values, with the opposite being true for higher values of input levels. • c =  = 1  identity function 97 Wolberg: Image Processing Course NotesGamma Correction • Cathode ray tube (CRT) devices γ (in) have an intensity tovoltage response that is a power function, 1 1 with  varying from 1.8 to 2.5 0 1 0 1 • This darkens the picture. gamma correction • Gamma correction ↓ is done by 1/γ γ (in) preprocessing the image before inputting it to the monitor. 1 1 1 0 0 1 98 Wolberg: Image Processing Course NotesExample: MRI (a) Dark MRI. Expand graylevel range for contrast manipulation  1 (b)  = 0.6, c=1 (c)  = 0.4 (best result) (d)  = 0.3 (limit of acceptability) When  is reduced too much, the image begins to reduce contrast to the point where it starts to have a “washedout” look, especially in the background 99 Wolberg: Image Processing Course NotesExample: Aerial Image Washedout image. Shrink graylevel range  1 (b)  = 3.0 (suitable) (c)  = 4.0 (suitable) (d)  = 5.0 (High contrast; the image has areas that are too dark; some detail is lost) 100 Wolberg: Image Processing Course NotesPiecewiseLinear Transformation Functions •Advantage: The form of piecewise functions can be arbitrarily complex •Disadvantage: Their specification requires considerably more user input 101 Wolberg: Image Processing Course NotesContrast Stretching • Low contrast may be due to poor illumination, a lack of dynamic range in the imaging sensor, or even a wrong setting of a lens aperture during acquisition. • Applied contrast stretching: (r ,s ) = (r ,0) and (r ,s ) = (r ,L1) 1 1 min 2 2 max 102 Wolberg: Image Processing Course NotesGraylevel Slicing 103 Wolberg: Image Processing Course NotesBitplane slicing • Highlighting the contribution made to total image appearance by specific bits Bitplane 7 One 8bit byte • Suppose each pixel is (most significant) represented by 8 bits • Higherorder bits contain the majority of the visually significant data • Useful for analyzing the relative importance played by Bitplane 0 each bit of the image (least significant) 104 Wolberg: Image Processing Course NotesExample • The (binary) image for bit plane 7 can be obtained by processing the input image with a thresholding graylevel transformation. Map all levels between 0 and 127 to 0 Map all levels between 129 and 255 to 255 An 8bit fractal image 105 Wolberg: Image Processing Course Notes8Bit Planes Bitplane 7 Bitplane 6 Bit Bit Bit plane 5 plane 4 plane 3 Bit Bit Bit plane 2 plane 1 plane 0 106 Wolberg: Image Processing Course NotesHardware LUTs • All point operations can be implemented by LUTs. • Hardware LUTs operate on the data as it is being displayed. • It’s an efficient means of applying transformations because changing display characteristics only requires loading a new table and not the entire image. • For a 1024x1024 8bit image, this translates to 256 entries instead of one million. • LUTs do not alter the contents of original image (nondestructive). Refresh memory Lookup table Display screen For display 0 0 0 1 20 20 20 20 40 V (i,j) V (i,j) 2 2 1 1 in 40 out 100 100 40 40 2 1 1 1 100 100 40 40 40 2 2 3 3 100 100 100 100 100 107 Wolberg: Image Processing Course NotesHistogram • A histogram of a digital image with gray levels in the range 0, L1 is a discrete function h(r ) = n k k th r : the k gray level k n : the number of pixels in the image having gray level r k k • The sum of all histogram entries is equal to the total number of pixels in the image. h(r) r 108 Wolberg: Image Processing Course NotesHistogram Example 5x5 image Plot of the Histogram Graylevel Count 6 0 2 5 2 3 4 4 6 1 5 4 2 3 1 2 4 5 6 3 4 3 1 1 5 6 6 4 5 2 5 2 0 1 3 3 4 1 6 4 0 1 2 3 4 Total 25 0 Pixel value Histogram evaluation: for(i=0; iMXGRAY; i++) Hi = 0; for(i=0; itotal; i++) Hini++; 109 Wolberg: Image Processing Course Notes CountNormalized Histogram • Divide each histogram entry at gray level r by k the total number of pixels in the image, n p( r ) = n / n k k • p( r ) gives an estimate of the probability of k occurrence of gray level r k • The sum of all components of a normalized histogram is equal to 1. 110 Wolberg: Image Processing Course NotesHistogram Processing •Basic for numerous spatial domain processing techniques. •Used effectively for image enhancement: Histogram stretching Histogram equalization Histogram matching •Information inherent in histograms also is useful in image compression and segmentation. 111 Wolberg: Image Processing Course NotesExample: Dark/Bright Images Dark image Components of histogram are concentrated on the low side of the gray scale. Bright image Components of histogram are concentrated on the high side of the gray scale. 112 Wolberg: Image Processing Course NotesExample: Low/High Contrast Images Lowcontrast image histogram is narrow and centered toward the middle of the gray scale Highcontrast image histogram covers broad range of the gray scale and the distribution of pixels is not too far from uniform, with very few vertical lines being much higher than the others 113 Wolberg: Image Processing Course NotesHistogram Stretching h(f) h(g) f g MIN MAX 255 0 3) Rescale to 0,255 range 1) Slide histogram down to 0 255( f MIN) g 2) Normalize histogram to 0,1 range MAX MIN 114 Wolberg: Image Processing Course NotesExample (1) Wide dynamic range permits for only a small improvement after histogram stretching 11 207 Image appears virtually identical to original 0 255 115 Wolberg: Image Processing Course NotesExample (2) • Improve effectiveness of histogram stretching by clipping intensities first Flat histogram: every graylevel is equally present in image 11 0 255 128 116 Wolberg: Image Processing Course NotesHistogram Equalization (1) • Produce image with flat histogram • All graylevels are equally likely • Appropriate for images with wide range of graylevels • Inappropriate for images with few graylevels (see below) 117 Wolberg: Image Processing Course NotesHistogram Equalization (2) Objective: we want a uniform histogram. Rationale: maximize image entropy. total h (v ) constant 1 out MXGRAY  h avg c (v ) (v1)h 1 out out avg This is a special case of histogram matching. Perfectly flat histogram: Hi = total/MXGRAY for 0 i MXGRAY. If Hv = k h then v must be mapped onto k different levels, avg from v to v . This is a oneto many mapping. 1 k 118 Wolberg: Image Processing Course NotesHistogram Equalization Mappings Rule 1: Always map v onto (v +v )/2. (This does not 1 k result in a flat histogram, but one where brightness levels are spaced apart). Rule 2: Assign at random one of the levels in v ,v . 1 k This can result in a loss of contrast if the original histogram had two distinct peaks that were far apart (i.e., an image of text). Rule 3: Examine neighborhood of pixel, and assign it a level from v ,v which is closest to neighborhood 1 k average. This can result in bluriness; more complex. Rule (1) creates a lookup table beforehand. Rules (2) and (3) are runtime operations. 119 Wolberg: Image Processing Course NotesExample (1) before after Histogram equalization 120 Wolberg: Image Processing Course NotesExample (2) before after Histogram equalization The quality is not improved much because the original image already has a wide graylevel scale 121 Wolberg: Image Processing Course NotesImplementation (1) No. of pixels 6 2 3 3 2 5 4 4 2 4 3 3 3 2 3 5 2 2 4 2 4 1 Gray level 4x4 image 0 1 2 3 4 5 6 7 8 9 Gray scale = 0,9 histogram 122 Wolberg: Image Processing Course NotesImplementation (2) Gray 0 1 2 3 4 5 6 7 8 9 Level(j) No. of 0 0 6 5 4 1 0 0 0 0 pixels k n  j 0 0 6 11 15 16 16 16 16 16 j0 k n j 6/16 11/16 15/16 16/16 16/16 16/16 16/16 16/16 0 0 s  n j0 3.3 6.1 8.4 0 0 9 9 9 9 9 s x 9 368 123 Wolberg: Image Processing Course NotesImplementation (3) No. of pixels 6 3 6 6 3 5 4 8 3 8 6 3 6 3 6 9 2 3 8 3 8 1 Output image 0 1 2 3 4 5 6 7 8 9 Gray scale = 0,9 Gray level Histogram equalization 124 Wolberg: Image Processing Course NotesNote (1) • Histogram equalization distributes the graylevels to reach maximum gray (white) because the cumulative distribution function equals 1 when 0  r  L1 k n  j • If is slightly different among consecutive k , those j0 graylevels will be mapped to (nearly) identical values as we have to produce an integer grayvalue as output • Thus, the discrete transformation function cannot guarantee a onetoone mapping 125 Wolberg: Image Processing Course NotesNote (2) • The implementation described above is widely interpreted as histogram equalization. • It is readily implemented with a LUT. • It does not produce a strictly flat histogram. • There is a more accurate solution. However, it may require a onetomany mapping that cannot be implemented with a LUT. 126 Wolberg: Image Processing Course NotesBetter Implementation (1) void histeq(imageP I1, imageP I2) int i, R; int leftMXGRAY, widthMXGRAY; uchar in, out; long total, Hsum, Havg, histoMXGRAY; / total number of pixels in image / total = (long) I1width I1height; / init I2 dimensions and buffer / I2width = I1width; I2height = I1height; I2image = (uchar ) malloc(total); / init input and output pointers / in = I1image; / input image buffer / out = I2image; / output image buffer / / compute histogram / for(i=0; iMXGRAY; i++) histoi = 0; / clear histogram / for(i=0; itotal; i++) histoini++; / eval histogram / R = 0; / right end of interval / Hsum = 0; / cumulative value for interval / Havg = total / MXGRAY; / interval value for uniform histogram / 127 Wolberg: Image Processing Course NotesBetter Implementation (2) / evaluate remapping of all input gray levels; Each input gray value maps to an interval of valid output values. The endpoints of the intervals are left and left+width. / for(i=0; iMXGRAY; i++) lefti = R; / left end of interval / Hsum += histoi; / cum. interval value / while(HsumHavg RMXGRAY1) / make interval wider / Hsum = Havg; / adjust Hsum / R++; / update right end / widthi = R lefti + 1; / width of interval / / visit all input pixels and remap intensities / for(i=0; itotal; i++) if(widthini == 1) outi = leftini; else / ini spills over into width possible values / / randomly pick from 0 to widthi / R = ((rand()0x7fff)widthini)15; / 0 = R width / outi = leftini + R; 128 Wolberg: Image Processing Course NotesNote • Histogram equalization has a disadvantage: it can generate only one type of output image. • With histogram specification we can specify the shape of the histogram that we wish the output image to have. • It doesn’t have to be a uniform histogram. • Histogram specification is a trialanderror process. • There are no rules for specifying histograms, and one must resort to analysis on a casebycase basis for any given enhancement task. 129 Wolberg: Image Processing Course NotesHistogram Matching C ( v ) 0 out C ( v’ ) 0 in C ( v’ ) 1 out h ( v ) h ( v ) 0 in 1 out v v in out ’ ’ v v in out v =T(v ) out in In the figure above, h() refers to the histogram, and c() refers to its cumulative histogram. Function c() is a monotonically increasing v function defined as: c(v) h(u)du  0 130 Wolberg: Image Processing Course NotesHistogram Matching Rule Let v =T(v ) If T() is a unique, monotonic function then out in v v out in h (u)du h (u)du 1 0  0 0 This can be restated in terms of the histogram matching rule: c (v ) c (v ) 1 out 0 in Where c (v ) = pixels  v , and c (v ) = pixels  v . 1 out out 0 in in. This requires that 1 v c (c (v )) out 1 0 in which is the basic equation for histogram matching techniques. 131 Wolberg: Image Processing Course NotesHistograms are Discrete • Impossible to match all histogram pairs because they are discrete. Refresh memory Lookup table Display screen For display 1 1 1 1 1 1 1 1 V (i,j) V (i,j) 1 1 2 2 in out 2 2 2 2 3 3 4 4 3 3 3 3 4 4 4 4 4 4 4 4 C ( v ) 1 out C ( v ) 0 in v v v out v in out in Continuous case Discrete case 132 Wolberg: Image Processing Course NotesProblems with Discrete Case • The set of input pixel values is a discrete set, and all the pixels of a given value are mapped to the same output value. For example, all six pixels of value one are mapped to the same value so it is impossible to have only four corresponding output pixels. 1 • No inverse for c in v = c (c (v )) because of 1 out 1 0 in discrete domain. Solution: choose v for which c (v ) out 1 out is closest to c (v ). 0 in • v v such that c (v ) c (v ) is a minimum in out 1 out 0 in 133 Wolberg: Image Processing Course NotesHistogram Matching Example (1) Input Histogram Input image Histogram match Output image Target Histogram 134 Wolberg: Image Processing Course NotesHistogram Matching Example (2) 135 Wolberg: Image Processing Course NotesImplementation (1) int histogramMatch(imageP I1, imageP histo, imageP I2) int i, p, R; int leftMXGRAY, rightMXGRAY; int total, Hsum, Havg, h1MXGRAY, h2; unsigned char in, out; double scale; / total number of pixels in image / total = (long) I1height I1width; / init I2 dimensions and buffer / I2width = I1width; I2height = I1height; I2image = (unsigned char ) malloc(total); in = I1image; / input image buffer / out = I2image; / output image buffer / for(i=0; iMXGRAY; i++) h1i = 0; / clear histogram / for(i=0; itotal; i++) h1ini++; / eval histogram / 136 Wolberg: Image Processing Course NotesImplementation (2) / target histogram / h2 = (int ) histoimage; / normalize h2 to conform with dimensions of I1 / for(i=Havg=0; iMXGRAY; i++) Havg += h2i; scale = (double) total / Havg; if(scale = 1) for(i=0; iMXGRAY; i++) h2i = scale; R = 0; Hsum = 0; / evaluate remapping of all input gray levels; Each input gray value maps to an interval of valid output values. The endpoints of the intervals are left and right / for(i=0; iMXGRAY; i++) lefti = R; / left end of interval / Hsum += h1i; / cumulative value for interval / while(Hsumh2R RMXGRAY1) / compute width of interval / Hsum = h2R; / adjust Hsum as interval widens / R++; / update / righti = R; / init right end of interval / 137 Wolberg: Image Processing Course NotesImplementation (3) / clear h1 and reuse it below / for(i=0; iMXGRAY; i++) h1i = 0; / visit all input pixels / for(i=0; itotal; i++) p = leftini; if(h1p h2p) / mapping satisfies h2 / outi = p; else outi = p = leftini = MIN(p+1, rightini); h1p++; 138 Wolberg: Image Processing Course NotesLocal Pixel Value Mappings • Histogram processing methods are global, in the sense that pixels are modified by a transformation function based on the graylevel content of an entire image. • We sometimes need to enhance details over small areas in an image, which is called a local enhancement. • Solution: apply transformation functions based on graylevel distribution within pixel neighborhood. 139 Wolberg: Image Processing Course NotesGeneral Procedure • Define a square or rectangular neighborhood. • Move the center of this area from pixel to pixel. • At each location, the histogram of the points in the neighborhood is computed and histogram equalization, histogram matching, or other graylevel mapping is performed. • Exploit easy histogram update since only one new row or column of neighborhood changes during pixelto pixel translation. • Another approach used to reduce computation is to utilize nonoverlapping regions, but this usually produces an undesirable checkerboard effect. 140 Wolberg: Image Processing Course NotesExample: Local Enhancement a) Original image (slightly blurred to reduce noise) b) global histogram equalization enhances noise slightly increases contrast but the structural details are unchanged c) local histogram equalization using 7x7 neighborhood reveals the small squares inside of the larger ones in the original image. 141 Wolberg: Image Processing Course NotesDefinitions (1) 1 mean (x, y) f (i, j)  n i, j 1 2 standard deviation  (x, y) ( f (i, j)(x, y))  n i, j • Let p(r ) denote the normalized histogram entry for grayvalue r i i for 0 ≤ i L where L is the number of graylevels. • It is an estimate of the probability of occurrence of graylevel r . i • Mean m can be rewritten as L1 m r p(r )  i i i0 142 Wolberg: Image Processing Course NotesDefinitions (2) • The nth moment of r about its mean is defined as L1 n  (r) (r m) p(r )  n i i i0 • It follows that: th 0 moment  (r) 1 0 1st moment  (r) 0 1 L1 2 nd  (r) (r m) p(r ) 2 moment 2 i i i0 2  (r) • The second moment is known as variance • The standard deviation is the square root of the variance. • The mean and standard deviation are measures of average grayvalue and average contrast, respectively. 143 Wolberg: Image Processing Course NotesExample: Statistical Differencing • Produces the same contrast throughout the image. • Stretch f(x, y) away from or towards the local mean to achieve a balanced local standard deviation throughout the image. • is the desired standard deviation and it controls the amount of stretch. 0 • The local mean can also be adjusted:  0 g(x, y)m (1)(x, y) ( f (x, y)(x, y)) 0  (x, y) • m is the mean to force locally and α controls the degree to which it is forced. 0 • To avoid problems when (x, y) = 0,  0 g(x, y)m (1)(x, y) ( f (x, y)(x, y)) 0   (x, y) 0 • Speedups can be achieved by dividing the image into blocks (tiles), exactly computing the mean and standard deviation at the center of each block, and then linearly interpolating between blocks in order to compute an approximation at any arbitrary position. In addition, the mean and standard deviation can be computed incrementally. 144 Wolberg: Image Processing Course NotesExample: Local Statistics (1) The filament in the center is clear. There is another filament on the right side that is darker and hard to see. Goal: enhance dark areas while leaving the light areas unchanged. 145 Wolberg: Image Processing Course NotesExample: Local Statistics (2) Solution: Identify candidate pixels to be dark pixels with low contrast. Dark: local mean k global mean, where 0 k 1. 0 0 Low contrast: k global variance local variance k global variance, 1 2 where k k . 1 2 Multiply identified pixels by constant E1. Leave other pixels alone. 146 Wolberg: Image Processing Course NotesExample: Local Statistics (3) Results for E=4, k =0.4, k =0.02, k =0.4. 3x3 neighborhoods used. 0 1 2 147 Wolberg: Image Processing Course NotesArithmetic/Logic Operations Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we describe arithmetic and logic operations commonly used in image processing. • Arithmetic ops: Addition, subtraction, multiplication, division Hybrid: crossdissolves • Logic ops: AND, OR, XOR, BIC, … 149 Wolberg: Image Processing Course NotesArithmetic/Logic Operations • Arithmetic/Logic operations are performed on a pixel bypixel basis between two images. • Logic NOT operation performs only on a single image. It is equivalent to a negative transformation. • Logic operations treat pixels as binary numbers: 158 235 = 10011110 11101011 = 10001010 • Use of LUTs requires 16bit rather than 8bit indices: Concatenate two 8bit input pixels to form a 16bit index into a 64Kentry LUT. Not commonly done. 150 Wolberg: Image Processing Course NotesAddition / Subtraction in1 +, , ,  out in2 Addition: for(i=0; itotal; i++) outi = MIN(((int)in1i+in2i), 255); Avoid overflow: clip result Subtraction: for(i=0; itotal; i++) outi = MAX(((int)in1iin2i), 0); Avoid underflow: clip result 151 Wolberg: Image Processing Course NotesOverflow / Underflow • Default datatype for pixel is unsigned char. • It is 1 byte that accounts for nonnegative range 0,255. • Addition of two such quantities may exceed 255 (overflow). • This will cause wraparound effect: 254: 11111110 255: 11111111 256: 100000000 257: 100000001 • Notice that loworder byte reverts to 0, 1, … when we exceed 255. • Clipping is performed to prevent wraparound. • Same comments apply to underflow (result 0). 152 Wolberg: Image Processing Course NotesImplementation Issues • The values of a subtraction operation may lie between 255 and 255. Addition: 0,510. • Clipping prevents over/underflow. • Alternative: scale results in one of two ways: 1. Add 255 to every pixel and then divide by 2. • Values may not cover full 0,255 range • Requires short intermediate image • Fast and simple to implement 2. Add negative of min difference (shift min to 0). Then, multiply all pixels by 255/(max difference) to scale range to 0,255 interval. • Full utilization of 0,255 range • Requires short intermediate image • More complex and difficult to implement 153 Wolberg: Image Processing Course NotesExample of Subtraction Operation 154 Wolberg: Image Processing Course NotesExample: Mask Mode Radiography • h(x,y) is the mask, an Xray image of a region of a patient’s body captured by an intensified TV camera (instead of traditional Xray film) located opposite an X ray source • f(x,y) is an Xray image taken after injection a contrast medium mask image h(x,y) image f(x,y) taken after into the patient’s bloodstream injection of a contrast medium (iodine) into the • images are captured at TV rates, bloodstream, with mask so the doctor can see how the subtracted out. Note: medium propagates through the • the background is dark because it doesn’t various arteries in an animation change much in both images. of f(x,y)h(x,y). • the difference area is bright because it has a big change 155 Wolberg: Image Processing Course NotesArithmetic Operations: CrossDissolve • Linearly interpolate between two images. • Used to perform a fade from one image to another. • Morphing can improve upon the results shown below. for(i=0; itotal; i++) outi = in1if + in2i(1f); f in 1 1 in 2 0 time in1 in2 156 Wolberg: Image Processing Course NotesMasking • Used for selecting subimages. • Also referred to as region of interest (ROI) processing. • In enhancement, masking is used primarily to isolate an area for processing. • AND and OR operations are used for masking. 157 Wolberg: Image Processing Course NotesExample of AND/OR Operation 158 Wolberg: Image Processing Course NotesDigital Halftoning Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review digital halftoning techniques to convert grayscale images to bitmaps: Unordered (random) dithering Ordered dithering Patterning Error diffusion 160 Wolberg: Image Processing Course NotesBackground • An 8bit grayscale image allows 256 distinct gray levels. • Such images can be displayed on a computer monitor if the hardware supports the required number of intensity levels. • However, some output devices print or display images with much fewer gray levels. • In these cases, the grayscale images must be converted to binary images, where pixels are only black (0) or white (255). • Thresholding is a poor choice due to objectionable artifacts. • Strategy: sprinkle blackandwhite dots to simulate gray. • Exploit spatial integration (averaging) performed by eye. 161 Wolberg: Image Processing Course NotesThresholding • The simplest way to convert from grayscale to binary. v out 0 255 v thr in 8 bpp (256 levels) 1 bpp (twolevel) Loss of information is unacceptable. 162 Wolberg: Image Processing Course NotesUnordered Dither (1) • Reduce quantization error by adding uniformly distributed white noise (dither signal) to the input image prior to quantization. • Dither hides objectional artifacts. • To each pixel of the image, add a random number in the range m, m, where m is MXGRAY/quantizationlevels. v out Uniform noise 0 255 v thr in 3 bpp (8 levels) 8 bpp (256 levels) 163 Wolberg: Image Processing Course NotesUnordered Dither (2) 1 bpp 2 bpp 3 bpp 4 bpp Quantization Dither/ Quantization 164 Wolberg: Image Processing Course NotesOrdered Dithering • Objective: expand the range of available intensities. • Simulates n bpp images with m bpp, where nm (usually m = 1). • Exploit eye’s spatial integration. Gray is due to average of black/white dot patterns. Each dot is a circle of black ink whose area is proportional to ( 1 – intensity). Graphics output devices approximate the variable circles of halftone reproductions. 1 0 2 3 4 • 2 x 2 pixel area of a bilevel display produces 5 intensity levels. 2 • n x n group of bilevel pixels produces n +1 intensity levels. • Tradeoff: spatial vs. intensity resolution. 165 Wolberg: Image Processing Course NotesDither Matrix (1) • Consider the following 2x2 and 3x3 dither matrices: 6 8 4  0 2   (2) (3) D D 1 0 3   3 1   5 2 7  • To display a pixel of intensity I, we turn on all pixels whose associated dither matrix values are less than I. • The recurrence relation given below generates larger dither matrices of dimension n x n, where n is a power of 2. (n / 2) (2) (n / 2) (n / 2) (2) (n / 2)  4D D U 4D D U (n) 00 01 D  (n / 2) (2) (n / 2) (n / 2) (2) (n / 2) 4D D U 4D D U 10 11  (n) where U is an n x n matrix of 1’s. 166 Wolberg: Image Processing Course NotesDither Matrix (2) • Example: a 4x4 dither matrix can be derived from the 2x2 matrix. 13 15 14 16 0 8 2 10  12  9 12 4 14 6 (4)  D  3 11 1 9  5 8 15 7 13 5  1 2 3 4 167 Wolberg: Image Processing Course NotesPatterning • Let the output image be larger than the input image. 2 • Quantize the input image to 0…n gray levels. • Threshold each pixel against all entries in the dither matrix. (4) Each pixel forms a 4x4 block of blackandwhite dots for a D matrix. An n x n input image becomes a 4n x 4n output image. • Multiple display pixels per input pixel. (n) • The dither matrix D is used as a spatiallyvarying threshold. ij • Large input areas of constant value are displayed exactly as before. 4n n 4n n 168 Wolberg: Image Processing Course NotesImplementation • Let the input and output images share the same size. 2 • First quantize the input image to 0…n gray levels. • Compare the dither matrix with the input image. for(y=0; yh; y++) // visit all input rows for(x=0; xw; x++) // visit all input cols i = x n; // dither matrix index j = y n; // dither matrix index (n) // threshold pixel using dither value D ij (n) outyw+x = (inyw+x D ) 255 : 0; ij 169 Wolberg: Image Processing Course NotesExamples 3 8 bpp (256 levels) 1 bpp (D ) 4 8 1 bpp (D ) 1 bpp (D ) 170 Wolberg: Image Processing Course NotesError Diffusion • An error is made every time a grayvalue is assigned to be black or white at the output. • Spread that error to its neighbors to compensate for over/undershoots in the output assignments If input pixel 130 is mapped to white (255) then its excessive brightness (255130) must be subtracted from neighbors to enforce a bias towards darker values to compensate for the excessive brightness. • Like ordered dithering, error diffusion permits the output image to share the same dimension as the input image. 171 Wolberg: Image Processing Course NotesFloydSteinberg Algorithm f(x, y) g(x, y) f (x, y) Threshold Distribute Error (w ) ij e(x, y) f (x, y) f (x, y) w e(x i, y j)"corrected intensity value"  ij i j  255 if f (x, y) MXGRAY / 2 g(x, y)  0 otherwise  e(x, y) f (x, y) g(x, y) w1  ij i j 172 Wolberg: Image Processing Course NotesError Diffusion Weights • Note that visual improvements are possible if lefttoright scanning among rows is replaced by serpentine scanning (zigzag). That is, scan odd rows from leftto right, and scan even rows from righttoleft. • Further improvements can be made by using larger neighborhoods. • The sum of the weights should equal 1 to avoid emphasizing or suppressing the spread of errors. x 7 / 48 5/ 48 x 8/ 42 4 / 42  x 7 /16 3/ 48 5/ 48 7 / 48 5/ 48 3/ 48 2 / 42 4 / 42 8/ 42 4 / 42 2 / 42 3/16 5 /16 1/16 1/ 48 3/ 48 5/ 48 3/ 48 1/ 48 1/ 42 2 / 42 4 / 42 2 / 42 1/ 42 FloydSteinberg JarvisJudiceNinke Stucki 173 Wolberg: Image Processing Course NotesExamples (1) FloydSteinberg JarvisJudiceNinke 174 Wolberg: Image Processing Course NotesExamples (2) FloydSteinberg JarvisJudiceNinke 175 Wolberg: Image Processing Course NotesExamples (3) FloydSteinberg JarvisJudiceNinke 176 Wolberg: Image Processing Course NotesImplementation thr = MXGRAY /2; // init threshold value for(y=0; yh; y++) // visit all input rows for(x=0; xw; x++) // visit all input cols out = (in thr) // threshold BLACK : WHITE; // note: use LUT e = in out; // eval error in 1 +=(e7/16.); // add error to E nbr inw1 +=(e3/16.); // add error to SW nbr in w +=(e5/16.); // add error to S nbr inw+1 +=(e1/16.); // add error to SE nbr in++; // advance input ptr out++; // advance output ptr 177 Wolberg: Image Processing Course NotesComments • Two potential problems complicate implementation: True for all errors can be deposited beyond image border neighborhood ops errors may force pixel grayvalues outside the 0,255 range Right border Right border  x 7 /16 x 7 / 48 5/ 48 Bottom border 3/16 5 /16 1/16 3/ 48 5/ 48 7 / 48 5/ 48 3/ 48 1/ 48 3/ 48 5/ 48 3/ 48 1/ 48 FloydSteinberg JarvisJudiceNinke 178 Wolberg: Image Processing Course NotesSolutions to Border Problem (1) • Perform if statement prior to every error deposit Drawback: inefficient / slow • Limit excursions of sliding weights to lie no closer than 1 pixel from image boundary (2 pixels for JJN weights). Drawback: output will be smaller than input • Pad image with extra rows and columns so that limited excursions will yield smaller image that conforms with original input dimensions. Padding serves as placeholder. Drawback: excessive memory needs for intermediate image input output padded input 179 Wolberg: Image Processing Course NotesSolutions to Border Problem (2) • Use of padding is further undermined by fact that 16 bit precision (short) is needed to accommodate pixel values outside 0, 255 range. • A better solution is suggested by fact that only two rows are active while processing a single scanline in the FloydSteinberg algorithm (3 for JJN). • Therefore, use a 2row (or 3row) circular buffer to handle the two (or three) current rows. • The circular buffer will have the necessary padding and 16bit precision. • This significantly reduces memory requirements. 180 Wolberg: Image Processing Course NotesCircular Buffer 0 0 0 1 1 1 2 2 3 3 2 4 4 1 5 5 2 3 4 3 input circular buffer output (snapshots) 181 Wolberg: Image Processing Course NotesNew Implementation thr = MXGRAY /2; // init threshold value copyRowToCircBuffer(0); // copy row 0 to circular buffer for(y=0; yh; y++) // visit all input rows copyRowToCircBuffer(y+1); // copy next row to circ buffer in1 = buf y 2 + 1; // circ buffer ptr; skip over pad in2 = buf(y+1)2 + 1; // circ buffer ptr; skip over pad for(x=0; xw; x++) // visit all input cols out = (in1 thr) BLACK : WHITE; // threshold e = in1 out; // eval error in1 1 +=(e7/16.); // add error to E nbr in21 +=(e3/16.); // add error to SW nbr in2 0 +=(e5/16.); // add error to S nbr in2 1 +=(e1/16.); // add error to SE nbr in1++; in2++ // advance circ buffer ptrs out++; // advance output ptr 182 Wolberg: Image Processing Course NotesNeighborhood Operations Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • This lecture describes various neighborhood operations: Blurring Edge detection Image sharpening Convolution 184 Wolberg: Image Processing Course NotesNeighborhood Operations • Output pixels are a function of several input pixels. • h(x,y) is defined to weigh the contributions of each input pixel to a particular output pixel. • g(x,y) = Tf(x,y); h(x,y) g(x,y)=Tf(x,y); h(x,y)=f(x,y)h(x,y) Filter h(x,y) Input: f(x,y) Output: g(x,y) 185 Wolberg: Image Processing Course NotesSpatial Filtering • h(x,y) is known as a filter kernel, filter mask, or window. • The values in a filter kernel are coefficients. • Kernels are usually of odd size: 3x3, 5x5, 7x7 • This permits them to be properly centered on a pixel Consider a horizontal crosssection of the kernel. Size of crosssection is odd since there are 2n+1 coefficients: n neighbors to the left + n neighbors to the right + center pixel h h h 1 2 3 h h h 4 5 6 h h h 7 8 9 186 Wolberg: Image Processing Course NotesSpatial Filtering Process • Slide filter kernel from pixel to pixel across an image. • Use raster order: lefttoright from the top to the bottom. • Let pixels have grayvalues f . i • The response of the filter at each (x,y) point is: 1D indexing R h f h f ... h f 1 1 2 2 mn mn mn 2D indexing  h f  i i ii Window centered at (4,3)  h ij g h f h f h f Kernel slides 43 1 32 2 33 3 34 across image  h f h f h f 4 42 5 43 6 44 In raster order  h f h f h f 7 52 8 53 9 54 187 Wolberg: Image Processing Course NotesLinear Filtering • Let f(x,y) be an image of size MxN. • Let h(i,j) be a filter kernel of size mxn. • Linear filtering is given by the expression: s t g(x, y) h(i, j) f (x i, y j)  is jt where s = (m1)/2 and t = (n1)/2 • For a complete filtered image this equation must be applied for x = 0, 1, 2, … , M1 and y = 0, 1, 2, … , N1. 188 Wolberg: Image Processing Course NotesSpatial Averaging • Used for blurring and for noise reduction • Blurring is used in preprocessing steps, such as removal of small details from an image prior to object extraction bridging of small gaps in lines or curves • Output is average of neighborhood pixels. • This reduces the “sharp” transitions in gray levels. • Sharp transitions include: random noise in the image edges of objects in the image • Smoothing reduces noise (good) and blurs edges (bad) 189 Wolberg: Image Processing Course Notes3x3 Smoothing Filters • The constant multiplier in front of each kernel is equal to the sum of the values of its coefficients. • This is required to compute an average. Weighted average Box filter The center is the most important and other pixels are inversely weighted as a function of their distance from the center of the mask. This reduces blurring in the smoothing process. 190 Wolberg: Image Processing Course NotesUnweighted/Weighted Averaging • Unweighted averaging (smoothing filter): 1 g(x, y) f (i, j)  m i, j • Weighted averaging: g(x, y) f (i, j)h(x i, y j)  i, j 7x7 unweighted averaging 7x7 Gaussian filter Original image 191 Wolberg: Image Processing Course NotesUnweighted Averaging • Unweighted averaging over a 5pixel neighborhood along a horizontal scanline can be done with the following statement: for(x=2; xw2; x++) outx=(inx2+inx1+inx+inx+1+inx+2)/5; • Each output pixel requires 5 pixel accesses, 4 adds, and 1 division. A simpler version (for unweighted averaging only) is: sum=in0+in1+in2+in3+in4; for(x=2; xw2; x++) outx = sum/5; sum+=(inx+3 – inx2); + Limited excursions reduce size of output 192 Wolberg: Image Processing Course NotesImage Averaging • Consider a noisy image g(x,y) formed by the addition of noise (x,y) to an original image f(x,y): g(x,y) = f(x,y) + (x,y) • If the noise has zero mean and is uncorrelated then we can compute the image formed by averaging K different noisy images: K 1 g(x, y) g (x, y)  i K i1 • The variance of the averaged image diminishes: 1 2 2  g (x, y)(x, y) K • Thus, as K increases the variability (noise) of the pixel at each location (x,y) decreases assuming that the images are all registered (aligned). 193 Wolberg: Image Processing Course NotesNoise Reduction (1) • Astronomy is an important application of image averaging. • Low light levels cause sensor noise to render single images virtually useless for analysis. 194 Wolberg: Image Processing Course NotesNoise Reduction (2) • Difference images and their histograms yield better appreciation of noise reduction. • Notice that the mean and standard deviation of the difference images decrease as K increases. 195 Wolberg: Image Processing Course NotesGeneral Form: Smoothing Mask • Filter of size mxn (where m and n are odd) s t h(i, j) f (x i, y j)  is jt g(x, y) s t h(i, j)  is jt summation of all coefficients of the mask Note that s = (m1)/2 and t = (n1)/2 196 Wolberg: Image Processing Course Notesa b Example c d e f • a) original image 500x500 pixel • b) f) results of smoothing with square averaging filter of size n = 3, 5, 9, 15 and 35, respectively. • Note: big mask is used to eliminate small objects from an image. the size of the mask establishes the relative size of the objects that will be blended with the background. 197 Wolberg: Image Processing Course NotesExample • Blur to get gross representation of objects. • Intensity of smaller objects blend with background. • Larger objects become bloblike and easy to detect. result of thresholding original image result after smoothing with 15x15 filter 198 Wolberg: Image Processing Course NotesUnsharp Masking • Smoothing affects transition regions where grayvalues vary. • Subtraction isolates these edge regions. • Adding edges back onto image causes edges to appear more pronounced, giving the effect of image sharpening. Blur + Sharpen Edge image image 199 Wolberg: Image Processing Course NotesOrderStatistics Filters • Nonlinear filters whose response is based on ordering (ranking) the pixels contained in the filter support. • Replace value of the center pixel with value determined by ranking result. • Order statistic filters applied to nxn neighborhoods: 2 median filter: R = medianz k = 1,2,…,n k 2 max filter: R = maxz k = 1,2,…, n k 2 min filter: R = minz k = 1,2,…, n k 200 Wolberg: Image Processing Course NotesMedian Filter • Sort all neighborhood pixels in increasing order. • Replace neighborhood center with the median. • The window shape does not need to be a square. • Special shapes can preserve line structures. • Useful in eliminating intensity spikes: salt pepper noise. 10 20 20 (10,15,20,20,20,20,25,25,200) 20 200 15 Median = 20 Replace 200 with 20 25 20 25 201 Wolberg: Image Processing Course NotesMedian Filter Properties • Excellent noise reduction • Forces noisy (distinct) pixels to conform to their neighbors. • Clusters of pixels that are light or dark with respect to their 2 neighbors, and whose area is less than n /2 (onehalf the filter area), are eliminated by an n x n median filter. • knearest neighbor is a variation that blends median filtering with blurring: Set output to average of k nearest entries around median 10 18 19 (10,15,18,19,20,20,25,25,200) k=1: replace 200 with (19+20+20)/3 20 200 15 k=2: replace 200 with (18+19+20+20+25)/5 25 20 25 k=3: replace 200 with (15+18+19+20+20+25+25)/7 k=4: replace 200 with (10+15+18+19+20+20+25+25+200)/9 202 Wolberg: Image Processing Course NotesExamples (1) Median filter output Additive salt pepper noise Blurring output 203 Wolberg: Image Processing Course NotesExamples (2) 204Derivative Operators • The response of a derivative operator is proportional to the degree of discontinuity of the image at the point at which the operator is applied. • Image differentiation enhances edges and other discontinuities (noise) deemphasizes area with slowly varying graylevel values. • Derivatives of a digital function are approximated by differences. 205FirstOrder Derivative • Must be zero in areas of constant grayvalues. • Must be nonzero at the onset of a grayvalue step or ramp. • Must be nonzero along ramps. f (x)  f (x1) f (x) x 206 Wolberg: Image Processing Course NotesSecondOrder Derivative • Must be zero in areas of constant grayvalues. • Must be nonzero at the onset of a grayvalue step or ramp. • Must be zero along ramps of constant slope. 2  f (x) f (x)f (x1) 2 x  f (x1) f (x1) 2 f (x) 207 Wolberg: Image Processing Course NotesExample 208 Wolberg: Image Processing Course NotesComparisons • 1storder derivatives: produce thicker edges strong response to graylevel steps • 2ndorder derivatives: strong response to fine detail (thin lines, isolated points) double response at step changes in graylevel 209 Wolberg: Image Processing Course NotesLaplacian Operator • Simplest isotropic derivative operator • Response independent of direction of the discontinuities. • Rotationinvariant: rotating the image and then applying the filter gives the same result as applying the filter to the image first and then rotating the result. • Since derivatives of any order are linear operations, the Laplacian is a linear operator. 2 2  f f 2  f 2 2 xy 210 Wolberg: Image Processing Course NotesDiscrete Form of Laplacian 2 2  f f 2  f 2 2 xy where 2  f  f (x1, y) f (x1, y) 2 f (x, y) 2 x 2  f  f (x, y1) f (x, y1) 2 f (x, y) 2 y 2  f f (x1, y) f (x1, y)  f (x, y1) f (x, y1) 4 f (x, y) 211 Wolberg: Image Processing Course NotesLaplacian Mask Isotropic result Isotropic result for rotations in for rotations in increments of 90° increments of 45° 212 Wolberg: Image Processing Course NotesAnother Derivation 1 1 1 Unweighted Average Smoothing Filter 1 1 1 1/9 1 1 1 0 0 0 Retain Original 0 1 0 0 0 0 1 1 1 Original – Average (negative of Laplacian Operator) 1 8 1 1/9 1 1 1 Summation of coefficients in masks equals 0. In constant areas: 0 Near edges: high values 213 Wolberg: Image Processing Course NotesEffect of Laplacian Operator • Since the Laplacian is a derivative operator it highlights graylevel discontinuities in an image it deemphasizes regions with slowly varying gray levels • The Laplacian tends to produce images that have grayish edge lines and other discontinuities all superimposed on a dark featureless background 214 Wolberg: Image Processing Course NotesExample 1 1 1 1 8 1 1 1 1 215 Wolberg: Image Processing Course NotesSimplification • Addition of image with Laplacian can be combined into one operator: g(x, y) f (x, y) f (x1, y) f (x1, y)  f (x, y1) f (x, y1) 4 f (x, y)  5 f (x, y) f (x1, y) f (x1, y)  f (x, y1) f (x, y1) 0 0 0 0 1 0 0 1 0 + 0 1 0 1 4 1 = 1 5 1 0 0 0 0 1 0 0 1 0 216 Wolberg: Image Processing Course NotesExample 217 Wolberg: Image Processing Course NotesGradient Operator (1) • The gradient is a vector of directional derivatives. f  f  x x f   f f y  y   • Although not strictly correct, the magnitude of the gradient vector is referred to as the gradient. • First derivatives are implemented using this magnitude f mag(f ) 1 2 2 approximation: 2  f f x y f f f 1 x y 2 2 2  ff     xy    218 Wolberg: Image Processing Course NotesGradient Operator (2) • The components of the gradient vector are linear operators, but the magnitude is not (square,square root). • The partial derivatives are not rotation invariant (isotropic), but the magnitude is. • The Laplacian operator yields a scalar: a single number indicating edge strength at point. • The gradient is actually a vector from which we can compute edge magnitude and direction. 2 2 f (i, j) f f or f (i, j) f f mag x y mag x y f y f (i, j) f (i1, j) f (i1, j) 1 x f (i, j) tan angle where f f (i, j) f (i, j1) f (i, j1) x y 219 Wolberg: Image Processing Course NotesSummary (1) Continuous Digital 1 1 1 1 f (x) v(i) 1 4 1 1 8 1 ' ' f (x) v (i) v(i) v(i1) 1 1 1 1 '' 2 '' ' ' f (x) f (x) v (i) v (i) v (i1) v(i) v(i1)v(i1) v(i 2)  v(i 2) 2v(i1) v(i) 1 2 1v(i 2) v(i1) v(i) 1 2 1v(i1) v(i) v(i1) The Laplacian is a scalar, giving only the magnitude about the change in pixel values at a point. The gradient gives both magnitude and direction. 220 Wolberg: Image Processing Course NotesSummary (2) Two dimensional : One dimensiona l : Sobel Operator : mask1 0 1 x 1 0 11 21  1   mask 2 0 2 mask 0 0 0  x y  mask 0 y   1 0 1 1 2 1   1  Prewitt Operator : 1 0 1111   mask1 0 1 mask 0 0 0 x y  1 0 1 1 1 1  Mag 2 2 in Gradient f (i, j) f f or f (i, j) f f mag x y mag x y Angle f y 1 f (i, j) tan angle f x Laplacian in Mag 221 Wolberg: Image Processing Course NotesExample (1) 222 Wolberg: Image Processing Course NotesExample (2) • Goal: sharpen image and bring out more skeletal detail. • Problem: narrow dynamic range and high noise content makes the image difficult to enhance. • Solution: 1. Apply Laplacian operator to highlight fine detail 2. Apply gradient operator to enhance prominent edges 3. Apply graylevel transformation to increase dynamic range 223 Wolberg: Image Processing Course Notes224 Wolberg: Image Processing Course Notes225 Wolberg: Image Processing Course NotesFiltering Theory Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • This lecture reviews filtering theory. Linearity and spatialinvariance (LSI) Impulse response Sifting integral Convolution 227 Wolberg: Image Processing Course NotesDefinitions • We use two criteria for filters: linearity and spatialinvariance g(x) g(x) f(x) f(x) Filter Linear : f (x)g(x) output is proportional to input f (x) f (x) g (x) g (x) superposition 1 2 1 2 or simply, f (x)f (x)g (x)g (x) 1 2 1 2 Space invariant(shift invariant) : f (x) g(x) 228 Wolberg: Image Processing Course NotesLSI Filter • Physically realizable filters (lenses) are rarely LSI filters. Optical systems impose a limit in their maximum response. Power can’t be negative, imposing a limit on the minimum response. Lens aberrations and finite image area prevent LSI property. • Nevertheless, close enough to approximate as LSI. 229 Wolberg: Image Processing Course NotesImpulse Response h(x) (x) Filter (x) Spike at x=0 (x)   ' ' lim f (x )dx1, x 0   (x)0    0, x 0   x x 0 0 (x) is the impulse function, or Dirac delta function which is a continuous function.  f (x ) f () (x)d  (x) samples f (x) at x 0 0 0 0   1, x 0  Kronecker delta function (discrete case: integer values of x)  (x)  0, x 0   Sifting integral f (x) f () (x)d   230 Wolberg: Image Processing Course NotesConvolution • Any input signal can be represented by an infinite sum of shifted and scaled impulses: (x) Input h(x) Filter x 0 0 Convolution: output Blurred version of input  continuous : g(x) f (x)h(x) f ()h(x)d   231 Wolberg: Image Processing Course NotesDiscrete Case • Output function is a scaled shifted version of impulse response. , : convolution operator h(x): convolution kernel, filter kernel • If h(x) = (x) then we have an ideal filter: output = input. • Usually h(x) extends over several neighbors. • Discrete convolution:  g(x) f (x) h(x) f ()h(x)   232 Wolberg: Image Processing Course NotesExample: Triangle Filter 1 1 x 0 x 1 h(x)  0 1 x  1 0 1 x Input samples (impulses) Impulse responses (additive) f 2 g f 1 g f (1 x) g f x 0 x 1 g 1 1 2 2 2 g 1 g g g ( f f )x f 1 2 2 1 1 0 1 x Addition of two adjacent impulse responses Straight line: y=mx+b • Linearity rule: scale h(x) according to f(x) and add g , g . 1 2 • Obtain f(x) for any x by sampling the reconstructed g(x). 233 Wolberg: Image Processing Course NotesConvolution Summation • g(x) is a continuous convolution summation to compute at particular values at x.  g(x) f (x) h(x) f ()h(x)   f 2 f 1 g(x) = f h(xx ) + f h(xx ) 1 1 2 2 If x =0 and x =1 then g(x)=f (1x)+f x 1 2 1 2 x x x 2 1 234 Wolberg: Image Processing Course NotesA Closer Look At the Convolution Integral f() h() f()h(x) f()h(x) 1 1/2  1 x x 1 1 1 x1 1 1 0 x 1 1 x 2 h(x) h() fh 1/2 1/2 1/2 x  1 1 2 1 235 Wolberg: Image Processing Course NotesMirrored Kernel • Why fold over (mirror) kernel h(x)  • Why not use g(x) f ()h( x)  f 2  f 1 By construction : f h(x x ) f h(x x ) 1 1 2 2 f ()h(x)  x x x 1 2 By centering h at x : f h(x x) f h(x x) 1 1 1 2 2 0 (which is wrong) Therefore flip h before centering at x. • We typically use symmetric kernels: h(x) = h(x) 1 1 0 1 x 236 Wolberg: Image Processing Course NotesImpulse Function • Impulse function (or Dirac delta function) is defined as   ' ' lim f (x )dx 1, x 0   (x)0    0, x 0   x 0 • It can be used to sample a continuous function f(x) at any point x , as follows: 0  f (x ) f () (x)d f (x ) (0) f (x ) 0 0 0 0    (x)  0 for  x 0 0 237 Wolberg: Image Processing Course NotesImpulse Response • When an impulse is applied to a filter, an altered impulse, (the impulse response) is generated at the (x) output. h(x) Filter x 0 • The first direct outcome of LSI is that the filter can be uniquely characterized by its impulse response. 238 Wolberg: Image Processing Course NotesSifting Integral • Any continuous input signal can be represented in the limit by an infinite sum of shifted and scaled impulses. • This is an outcome of the sifting integral:  f (x) f () (x)d   which uses signal f(x) to scale the collection of impulses: = 239 Wolberg: Image Processing Course NotesConvolution Integral (1) • The response g(x) of a digital filter to an arbitrary input signal f(x) is expressed in terms of the impulse response h(x) of the filter by means of convolution integral:  g(x) f (x)h(x) f ()h(x)d   where denotes the convolution operation, h(x) is used as the convolution (filter) kernel, and  is the dummy variable of integration. • Kernel h(x) is treated as a sliding window that is shifted across the entire input signal. • As h(x) makes its way across f(x), a sum of the pointwise products between the two functions is taken and assigned to output g(x). 240 Wolberg: Image Processing Course NotesConvolution Integral (2) • This process, known as convolution, is of fundamental importance to linear filtering theory. • It simply states that the output of an LSI filter will be a superposition of shifted and scaled impulse responses. • This is used to explain how a continuous image is blurred by a camera as it passes through the lens. In this context, h(x) is known as the point spread function (PSF), reflecting the limitation of the camera to accurately resolve a small point without somewhat blurring it. output Filter 241 Wolberg: Image Processing Course NotesEvaluation • We can arrive at g(x) in two ways: Graphical construction of a series of shifted/scaled impulse responses Computing the convolution summation at all points of interest • Graphical construction more closely follows the physical process as an input signal passes through a filter. • Convolution summation more closely follows the practical evaluation of g(x) at a finite number of desired points. instead of adding scaled and shifted unit impulse responses (responses of unit impulses), we center the unit impulse response at the point of interest. 242 Wolberg: Image Processing Course NotesGraphical Construction (x) Asymmetric unit impulse response 1 h(x) with an extent of 5 points .5 g(x ) Filter 0 2 x 0 2 0 1 2 1 1.5 A f =4 2 1 Blowup of f =3 4 1 B interval g(x ) C f =2 0 f =2 3 1 2 between D Filter 1.5 1 and 2 1 1 0 1 2 x 0 .5 h(x 2) 1 2 0 1 2 3 0 2 3 1 5 0 4 x 0 h(x 1) 0 g(x ) f h(x ) f h(x1) f h(x 2) f h(x 3) h(x 3) 0 0 1 0 2 0 3 0 4 0 h(x ) 0 g(x ) 2h(x ) 4h(x1) 2h(x 2) 3h(x 3) 0 0 0 0 0 x x 1 x 3 x 2 0 0 0 0 As we scan f from left to right we multiply samples with values taken from h i in right to left order. 243 Wolberg: Image Processing Course NotesConvolution Summation (1) • Instead of adding scaled/shifted unit impulse responses, center unit impulse response at point of interest: f =4 2 f =3 4 f =2 f =2 3 1 h(x ) 0  .5 g(x ) f ()h(x)d 0 0   0 1 x 2 3 0 • As  increases, we scan f() from left to right. However, h(x ) is scanned from right to left. 0 • Reason: points to the left of x had contributed to x 0 0 through the right side of h (see graphical construction). 244 Wolberg: Image Processing Course NotesConvolution Summation (2) • More straightforward: implement convolution if both the input and the kernel were scanned in the same direction. • This permits direct pointwise multiplication among them. • Thus, flip kernel before centering it on output position. f =4 2 f =3 4 f =2 3 f =2 1 .5 h 3 h 2 h 1 h 4 0 1 x 2 3 0 x x x x 1 2 3 4 245 Wolberg: Image Processing Course NotesConvolution Summation (3) 4 f h • will yield the value for g(x ), where h = h(x x ).  i i 0 i i 0 i1 • Rationale: multiplying the unit impulse response h with f i, h is being scaled. The distance between x and x 0 i accounts for the effect of a shifted response function on the current output. • For most impulse response functions, h will taper off with increasing distance from its center. • If the kernel is symmetric (h(x)=h(x)), it is not necessary to flip the kernel before centering. 246 Wolberg: Image Processing Course NotesDiscrete Convolution (1)  for integer values of  and arbitrary values of x g(x) f ()h(x)   • The kernel is often a discrete set of weights, i.e., a 3x3 filter kernel. • As long as kernel is shifted in pixel (integer) increments across the image, there is no alignment problem with underlying image. • However, for noninteger pixel increments the kernel values may have f no corresponding image values. 5 f 2 f 4 f f 1 3 h 0 h is aligned with f h h 1 1 h h 2 2 (h is centered at pixel f ) 3 h 0 h h 1 1 h h is not aligned with f h 2 2 It has been shifted by a ½pixel increment 247 Wolberg: Image Processing Course NotesDiscrete Convolution (2) There are two possible solutions to this problem: 1. Represent h in analytic form for evaluation anywhere 2 4x Ex: Let bellshape PSF be h(x)=2 . The appropriate weight to apply to f can be obtained by setting x to the difference between the i center of the bell (impulse response) and the position of the f . i 2. Supersample h so that a dense set of samples are used to represent h. The supersampled version will then be aligned with the image data, or at least it will have values that are nearby. f 5 f 2 f 4 f f If the kernel is known to be slid across the 1 3 image at fixed increments, then the kernel can be sampled at known positions. Supersampled kernel Ex: a 1/3 increment requires the kernel to Be sampled three times per unit interval. Closest kernel samples from supersampled set 248 Wolberg: Image Processing Course NotesFourier Transform Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • This lecture reviews Fourier transforms and processing in the frequency domain. Definitions Fourier series Fourier transform Fourier analysis and synthesis Discrete Fourier transform (DFT) Fast Fourier transform (FFT) 250 Wolberg: Image Processing Course NotesBackground (1) • Fourier proved that any periodic function can be expressed as the sum of sinusoids of different frequencies, each multiplied by a different coefficient. → Fourier series • Even aperiodic functions (whose area under the curve is finite) can be expressed as the integral of sinusoids multiplied by a weighting function. → Fourier transform • In a great leap of imagination, Fourier outlined these results in a memoir in 1807 and published them in La Theorie Analitique de la Chaleur (The Analytic theory of Heat) in 1822. The book was translated into English in 1878. 251 Wolberg: Image Processing Course NotesBackground (2) • The Fourier transform is more useful than the Fourier series in most practical problems since it handles signals of finite duration. • The Fourier transform takes us between the spatial and frequency domains. • It permits for a dual representation of a signal that is amenable for filtering and analysis. • Revolutionized the field of signal processing. 252 Wolberg: Image Processing Course NotesExample 253 Wolberg: Image Processing Course NotesUseful Analogy • A glass prism is a physical device that separates light into various color components, each depending on its wavelength (or frequency) content. • The Fourier transform is a mathematical prism that separates a function into its frequency components. 254 Wolberg: Image Processing Course NotesFourier Series for Periodic Functions 255 Wolberg: Image Processing Course NotesFourier Transform for Aperiodic Functions 256 Wolberg: Image Processing Course NotesFourier Analysis and Synthesis • Fourier analysis: determine amplitude phase shifts • Fourier synthesis: add scaled and shifted sinusoids together • Fourier transform pair:  i2ux Forward F.T. F(u) f (x)e dx    i2ux Inverse F.T. f (x) F(u)e dx   where i1, and i2ux e cos(2ux) i sin(2ux) complex exponential at freq. u Euler’s formula 257 Wolberg: Image Processing Course NotesFourier Coefficients • Fourier coefficients F(u) specify, for each frequency u, the amplitude and phase of each complex exponential. • F(u) is the frequency spectrum. • f(x) and F(u) are two equivalent representations of the same signal. F(u) R(u) iI(u) 2 2 F(u) R (u) I (u) magnitude spectrum; aka Fourier spectrum I(u) 1 (u) tan phase spectrum R(u) 2 P(u) F(u) 2 2  R (u) I (u) spectral density 258 Wolberg: Image Processing Course Notes1D Example A 0 x W A f (x)  0 x W  x W  1 i2ux ax ax F(u) f (x)e dx note : e dx e  a  W i2ux F(u) Ae dx  0 AW  A W A A i2uxi2uwi2uw  F(u) e e1 1 e 0 i2u i2u i2u ixix A e e iuwiuwiuw F(u) e e e note : sin x i2u 2i 1/w 3/w u 3/w 1/w 0 2/w 2/w A iuW F(u) sinuWe complex function u 1 A sinuW iuW F(u) sinuW e AW AW sin cuW sinc uuW sin(x) sin( x) where sinc(x) (some books have sinc(x) ) 0 4 2 2 4 x x 259 Wolberg: Image Processing Course Notes2D Example 260 Wolberg: Image Processing Course NotesFourier Series (1) For periodic signals, we have the Fourier series: n i2u x th 0 f (x) c(nu )e where c(nu ) is the n Fourier coefficient  0 0 n  x / 2 0 1 i2u x 0 c(nu ) f (x)e dx 0  x 0 x / 2 0 That is, the periodic signal contains all the frequencies that are harmonics of the fundamental frequency. 261 Wolberg: Image Processing Course NotesFourier Series (2) x / 2 W / 2 0 1 1 i2nu xi2nu x 0 0 c(nu ) f (x)e dx Ae dx 0  x x 0 0 x / 2W / 2 0 A inu Winu W 0 0 c(nu ) (e e ) 0  i2nu x 0 0 ixix A e e c(nu ) sin(nu W )  sin x ; u x 1 0 0 0 0 n 2i Au W 0 c(nu ) sin(nu W ) Au Wsinc(nu W ) 0 0 0 0 nu W 0 x W 0 Note that if  , then we have a square wave and 2 2 c(nu ) Au x sinc(nu x ) 0 0 0 0 0 A sinc(n) n1,3,  c(nu )  0 0 n 0,2,4,  262 Wolberg: Image Processing Course NotesFourier Series (3) • The Fourier transform is applied for aperiodic signals. • It is represented as an integral over a continuum of frequencies. • The Fourier Series is applied for periodic signals. • It is represented as a summation of frequency components that are integer multiples of some fundamental frequency. 263 Wolberg: Image Processing Course NotesExample Ex : Rectangular Pulse Train W  A x  2 f (x) in interval W/2, W/2  W  0 x  2 A W W time W/2 0 W/2 264 Wolberg: Image Processing Course NotesDiscrete Fourier Transform ux N1 i2 1 N F(u) f (x)e forward DFT  N x0 ux N1 i2 N f (x) F(u)e inverse DFT  u0 for 0 ≤ u ≤ N1 and 0 ≤ x ≤ N1 where N is the number of equispaced input samples. The 1/N factor can be in front of f(x) instead. 265 Wolberg: Image Processing Course NotesFourier Analysis Code • DFT maps N input samples of f into the N frequency terms in F. for(u=0; uN; u++) /compute spectrum over all freq. u / real = imag = 0; /reset real, imag component of F(u)/ for(x=0; xN; x++) / visit each input pixel / real += (fxcos(2PIux/N)); imag += (fxsin(2PIux/N)); / Note: if f is complex, then real += (frxcos()fixsin()); imag += (frxsin()+fixcos()); because (f +if )(g +ig )=(f g f g )+i(f g +f g ) r i r i r r i i i r r i / Fru = real / N; Fiu = imag / N; 266 Wolberg: Image Processing Course NotesFourier Synthesis Code for(x=0; xN; x++) / compute each output pixel / real = imag = 0; / reset real, imaginary component / for(u=0; uN; u++) c = cos(2PIux/N); s = sin(2PIux/N); real += (FrucFius); imag += (Frus+Fiuc); frx = real; / OR fx = sqrt(realreal + imagimag); fix = imag; 267 Wolberg: Image Processing Course NotesExample: Fourier Analysis (1) 1.5 1.5 1 1 0.5 0.5 0 0 0.5 0.5 1 1 1.5 1.5 0 2 4 6 8 10 0 2 4 6 8 10 t t 1.5 1.5 1 1 0.5 0.5 0 0 0.5 0.5 1 1 1.5 1.5 0 2 4 6 8 10 0 2 4 6 8 10 t t 268 Wolberg: Image Processing Course Notes square signal, sw(t) square signal, sw(t) square signal, sw(t) square signal, sw(t)Gibbs phenomenon Example: Fourier Analysis (2) 1.5 1.5 1 1 0.5 0.5 0 0 0.5 0.5 1 1 1.5 1.5 0 2 4 6 8 10 0 2 4 6 8 10 t t 1.5 1.5 1 1 0.5 0.5 0 0 0.5 0.5 1 1 1.5 1.5 0 2 4 6 8 10 0 2 4 6 8 10 t t 269 Wolberg: Image Processing Course Notes square signal, sw(t) square signal, sw(t) square signal, sw(t) square signal, sw(t)Summary 2.5 2 1.5 T 1 1 ikω t 0.5 Periodic c s(t)e dt k 0 FS Discrete 0 1 2 3 4 5 6 7 8 T 0 time, t (period T) 2.5 Continuous 2 i2πf t  1.5 S(f) s(t)e dt 1 Aperiodic Continuous  FT  0.5 0 0 2 4 6 8 10 12 time, t 2.5 2πk n N1 i 2 1 N 1.5 c sne  k 1 N Periodic Discrete DFS n0 0.5 0 0 1 2 3 4 5 6 7 8 (period T) time, t k  i2πf n Discrete S(f) sne  Continuous 2.5 DTFT n 2 1.5 Aperiodic 1 2πk n N1 0.5i 1 N 0 c sne 0 2 4 6 8 10 12 DFT k time, t Discrete k N n0 Note: i=1,  = 2/T, sn=s(t ), N = of samples n 270 Wolberg: Image Processing Course Notes2D Fourier Transform Continuous: i2 (uxvy)i2uxi2vy F f (x, y) F(u,v) f (x, y)e dx f (x, y)e e dx   1i2 (uxvy)i2uxi2vy F F(u,v) f (x, y) F(u,v)e dx F(u,v)e e dx   Separable: F(u, v) = F(u)F(v) Discrete: ux vy ux vy N1 M1 M1 N1 i2 ( )i2 ( )i2 ( ) 1 1 1 N M N M F(u,v) f (x, y)e dx f (x, y)e e  MN M N x0 y0 y0 x0  ux vy ux vy N1 M1 M1 N1 i2 ( )i2 ( )i2 ( ) N M N M f (x, y) F(u,v)e dx F(u,v)e e  x0 y0 y0 x0  271 Wolberg: Image Processing Course NotesSeparable Implementation N1 M1 1 1  j2vy / N j2ux / M F(u,v) e f (x, y)e  N M y0 x0 N1 1  j2vy / N transform each row  F(u, y)e  N y0 transform each column of intermediate result M1 1  j2ux / M where F(u, y) f (x, y)e  M x0 The 2D Fourier transform is computed in two passes: 1) Compute the transform along each row independently. 2) Compute the transform along each column of this intermediate result. 272 Wolberg: Image Processing Course NotesProperties • Edge orientations in image appear in spectrum, rotated by 90°. • 3 orientations are prominent: 45°, 45°, and nearly horizontal long white element. 273 Wolberg: Image Processing Course NotesMagnitude and Phase Spectrum 2DFourier Phase Mad.bw Magnitude transforms example Phase Mandrill.bw Magnitude 274 Wolberg: Image Processing Course NotesRole of Magnitude vs Phase (1) Rick Linda Pictures reconstructed using the Fourier phase of another picture MagLinda MagRick PhaseRick PhaseLinda 275 Wolberg: Image Processing Course NotesRole of Magnitude vs Phase (2) Magnitude Phase 276 Wolberg: Image Processing Course NotesNoise Removal 277 Wolberg: Image Processing Course NotesFast Fourier Transform (1) k nk W • The DFT was defined as: n ux N1 i2 1 N F(u) f (x)e 0 x N1  N x0 Rewrite: N nk N1 i2 1 freq. N F f e 0 n N1  n k N k0 N1 nk Let F f W  n k k0 i2  2 2 N where W e cos( ) i sin( ) N N r Also, Let N 2 (N is apower of 2) 278 Wolberg: Image Processing Course NotesFast Fourier Transform (2) k nk W nk • W can be thought of as a 2D n array, indexed by n and k. • It represents N equispaced values along a sinusoid at each of N frequencies. N freq. • For each frequency n, there are N multiplications (N samples in sine wave of freq. n). Since there are N frequencies, DFT: 2 O(N ) • With the FFT, we will derive an O(N logN) process. 279 Wolberg: Image Processing Course NotesComputational Advantage 280 Wolberg: Image Processing Course NotesDanielson–Lanczos Lemma (1) 1942: nk N1 i2 N F f e n k k0 N N 11 n(2k ) n(2k1) 2 2 i2i2 N N F f e f e n 2k 2k1 k0 k0 Even Numbered Terms Odd Numbered Terms f , f , f , f , f , f , 0 2 4 1 3 5 281 Wolberg: Image Processing Course NotesDanielson–Lanczos Lemma (2) N N 11 i2nki2nk 2 2 N 2 n N 2 F f eW f e n 2k 2k1 k0 k0 i2 N W e e n o F FW F n n n th n component of F.T. of length N/2 formed from the even components of f th n component of F.T. of length N/2 formed from the odd components of f DivideandConquer solution: Solving a problem (F ) is reduced to 2 n smaller ones. e o Potential Problem: n in F and F is still made to vary from 0 to N1. Since n n each subproblem is no smaller than original, it appears wasteful. Solution: Exploit symmetries to reduce computational complexity. 282 Wolberg: Image Processing Course NotesDanielson–Lanczos Lemma (3) Given : a DFT of length N, F F nN n i2 (nN )ki2nki2Nk N1 N1 N N N Proof : F f e f e e nN k k k0 k0 i2nki2nk N1 N1 N N  f e (cos2k i sin 2k) f e F  k k n k0 k0 N n  2 N 2 N  2 W cos n i sin n   N 2 N 2    2n 2n   cos i sin  N N   2n 2n  cos i sin  N N  n W 283 Wolberg: Image Processing Course NotesMain Points of FFT nk N1 i2 N F f e n k k0 e n o F F W F n n n N N length N length length 2 2 but 0 n N = e n o F F + W F e n n n F n n ±W = e n o F F W F o n n n F n N N/2 284 Wolberg: Image Processing Course NotesFFT Example (1) • Input: 10, 15, 20, 25, 5, 30, 8, 4 10 15 20 25 5 30 8 4 e o 10 20 5 8 15 25 30 4 e o e o 10 5 20 8 15 30 25 4 e o e o e o e o 10 5 20 8 15 30 25 4 eee eeo eoe eoo oee oeo ooe ooo 000 001 010 011 100 101 110 111 285 Wolberg: Image Processing Course NotesFFT Example (2) 10 5 20 8 15 30 25 4 0 0 0 0 = e n o = e n o F F + W F F F W F n n n n n n 15 5 28 12 45 15 29 21 1i2 / 4 0 1 0 1 W ei 5 5+ 15 15 43 13 74 16 12i 12i 21i +21i 1i2 / 8 W e B 0 1 2 3 2i4 / 8 W e C X=(512i)+ 3i6 / 8 (1521i)B W e D 13 13 Y X (see next page for X Y 31 117 Y=(5+12i)+ 16i +16i weights B, C, D) (15+21i)D 286 Wolberg: Image Processing Course NotesWeights i2πux/N • DFT is a convolution with kernel values e • These values are derived from a unit circle. i G 2 2 2 2  i  i F H 2 2 2 2 E A 1 1 B D 2 2 2 2  i  i C 2 2 2 2 i 287 Wolberg: Image Processing Course NotesDFT Example (1) i 2 2 2 2 G  i  i 2 2 • Input: 10, 15, 20, 25, 5, 30, 8, 4 2 2 F H nk E A 1 1 N1 i2 N F f e  n k 2 2 2 2 B D  i i k0 2 2 C 2 2 i F 1015 20 25 5 30 8 4 117 0 i2 (1)(0) / 8i2 (1)(1) / 8i2 (1)(2) / 8i2 (1)(6) / 8i2 (1)(7) / 8 F 10e15e 20e ... 8e 4e 1  10A15B 20C 25D 5E 30F 8G 4H i2 (2)(0) / 8i2 (2)(1) / 8i2 (2)(2) / 8i2 (2)(6) / 8i2 (2)(7) / 8 F 10e15e 20e ... 8e 4e 2  10A15C 20E 25G 5A 30C 8E 4G1316i i2 (3)(0) / 8i2 (3)(1) / 8i2 (3)(2) / 8i2 (3)(6) / 8i2 (3)(7) / 8 F 10e15e 20e ... 8e 4e 3  10A15D 20G 25B 5E 30H 8C 4F 288 Wolberg: Image Processing Course Notesi 2 2 2 2 G  i  i DFT Example (2) 2 2 2 2 F H A E 1 1 2 2 2 2 B D  i  i 2 2 C 2 2 i i2 (4)(0) / 8i2 (4)(1) / 8i2 (4)(2) / 8i2 (4)(6) / 8i2 (4)(7) / 8 F 10e15e 20e ... 8e 4e 4  10A15E 20A 25E 5A 30E 8A 4E31 i2 (5)(0) / 8i2 (5)(1) / 8i2 (5)(2) / 8i2 (5)(6) / 8i2 (5)(7) / 8 F 10e15e 20e ... 8e 4e 5  10A15F 20C 25H 5E 30B 8G 4D i2 (6)(0) / 8i2 (6)(1) / 8i2 (6)(2) / 8i2 (6)(6) / 8i2 (6)(7) / 8 F 10e15e 20e ... 8e 4e 6  10A15G 20E 25C 5A 30G 8E 4C1316i i2 (7)(0) / 8i2 (7)(1) / 8i2 (7)(2) / 8i2 (7)(6) / 8i2 (7)(7) / 8 F 10e15e 20e ... 8e 4e 7  10A15H 20G 25F 5E 30D 8C 4B 289 Wolberg: Image Processing Course Notes290 Wolberg: Image Processing Course Notes291 Wolberg: Image Processing Course Notes292 Wolberg: Image Processing Course Notes293 Wolberg: Image Processing Course NotesFiltering in the Frequency Domain Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • This lecture reviews frequency domain filtering. Convolution theorem Frequency bands Lowpass filter • Ideal, Butterworth, Gaussian Highpass filter • Ideal, Butterworth, Gaussian Notch filter Homomorphic filtering 295 Wolberg: Image Processing Course NotesBasic Steps 296 Wolberg: Image Processing Course NotesBasics x+y 1. Multiply input image by (1) to center the transform to u = M/2 and v = N/2 (if M and N are even numbers, then shifted coordinates will be integers) 2. Compute F(u,v), the DFT of the image from (1) 3. Multiply F(u,v) by a filter function H(u,v) 4. Compute the inverse DFT of the result in (3) 5. Obtain the real part of the result in (4) x+y 6. Multiply the result in (5) by (1) to cancel the multiplication of the input image. 297 Wolberg: Image Processing Course NotesConvolution Theorem “Multiply F(u,v) by a filter function H(u,v)” This step exploits the convolution theorem: Convolution in one domain is equivalent to multiplication in the other domain. f(x)g(x) ↔ F(u) G(u) f(x) g(x) ↔ F(u)G(u) 298 Wolberg: Image Processing Course NotesPeriodicity F(u,v) = F(u+M, v) = F(u,v+N) = F(u+M,v+N) 299 Wolberg: Image Processing Course NotesComments • The magnitudes from (M/2)+1 to M1 are reflections of the values in the half period to the left of the origin. • DFT is formulated for values of u in the interval 0,M1. • This yields two backtoback half periods in this interval. • For one full period, move origin to u=M/2 by multiplying x f(x) by (1) . 300 Wolberg: Image Processing Course NotesSignificance of Periodicity x ranges from 0 to 799 period too short for h() to slide past f 301 Wolberg: Image Processing Course NotesPadding (1) • Corrupted results were due to inadequate period. • Solution: add padding to increase period. 302 Wolberg: Image Processing Course NotesPadding (2) 303 Wolberg: Image Processing Course NotesExample 304 Wolberg: Image Processing Course NotesFrequency Bands • Low frequencies: general graylevel appearance over smooth areas: slowly varying grayscales. • High frequencies: responsible for detail, such as edges and noise. • A filter that attenuates high frequencies while “passing” low frequencies is a lowpass filter. • A filter that attenuates low frequencies while “passing” high frequencies is a highpass filter. • Lowpass filters blur images. • Highpass filters highlight edges. 305 Wolberg: Image Processing Course NotesLowpass and Highpass Filters Lowpass filter Highpass filter 306 Wolberg: Image Processing Course NotesImproved Highpass Output • Add constant to filter so that it will not eliminate F(0,0) 307 Wolberg: Image Processing Course NotesCorresponding Filters in the Spatial and Frequency Domains 308 Wolberg: Image Processing Course NotesNotch Filter • F(0,0) represents the average value. • If F(0,0) = 0, then average will be 0. • Achieved by applying a notch filter: 0 if (u, v)  (M/2, N/2 )  H (u,v)  1 otherwise  • In reality the average of the displayed image can’t be zero as it needs to have negative gray levels. The output image needs to scale the graylevel. • Filter has notch (hole) at origin. 309 Wolberg: Image Processing Course NotesIdeal Lowpass Filter cutoff frequency 310 Wolberg: Image Processing Course NotesImage Power Circles 311 Wolberg: Image Processing Course NotesVarying ILPF Cutoff Frequencies • Most sharp detail in this image is contained in the 8 power removed by filter. • Ringing behavior is a characteristic of ideal filters. • Little edge info contained in upper 0.5 of spectrum power in this case. 312 Wolberg: Image Processing Course NotesILPF Ringing • Ringing in spatial filter has two characteristics: dominant component at origin and concentric circular components. • Center component responsible for blurring. • Concentric components responsible for ringing. • Radius of center comp. and circles/unit distance are inversely proportional to cutoff frequency. 313 Wolberg: Image Processing Course NotesButterworth Lowpass Filter (BLPF) 1 H (u,v) n is filter order 2n 1 D(u,v) / D 0 314 Wolberg: Image Processing Course NotesVarying BLPF Cutoff Frequencies • Unlike the ILPF, BLPF does not have a sharp discontinuity at the cutoff frequency. • At D , H(u,v) is down 50 0 from max value of 1. • Smooth transition in blurring as cutoff frequency increases. • No ringing due to filter’s smooth transition between low and high frequencies. 315 Wolberg: Image Processing Course NotesSpatial Representation of BLPFs No ringing Imperceptible ringing Significant ringing 316 Wolberg: Image Processing Course NotesGaussian Lowpass Filter (GLPF) 2 2 σ = D = cutoff freq D (u,v) / 2D 0 0 H(u,v) e 317 Wolberg: Image Processing Course NotesVarying GLPF Cutoff Frequencies • Smooth transition in blurring as cutoff frequency increases. • GLPF did not achieve as much smoothing as BLPF of order 2 for same cutoff freq. • GLPF profile is not as tight as profile of the BLPF of order 2. • No ringing. • BLPF is more suitable choice if tight control is needed around cutoff frequency. Only drawback: possible ringing. 318 Wolberg: Image Processing Course NotesExample (1) 319 Wolberg: Image Processing Course NotesExample (2) 320 Wolberg: Image Processing Course NotesExample (3) • Crude but simple way of removing prominent scan lines. • Blurs out detail while leaving large features recognizable. • Averages out features smaller than the ones of interest. 321 Wolberg: Image Processing Course NotesHighpass Filters: 1H LP Ideal highpass filter 0 if D(u, v) D  0 H (u,v)  1 if D(u, v) D 0  Butterworth highpass filter 1 H(u,v) 2n 1D D(u,v) 0 Gaussian highpass filter 2 2 D (u,v)/ 2D 0 H(u,v)1 e 322 Wolberg: Image Processing Course NotesSpatial Representation of Filters Ideal Butterworth Gaussian 323 Wolberg: Image Processing Course NotesVarying IHPF Cutoff Frequencies 324 Wolberg: Image Processing Course NotesVarying BHPF Cutoff Frequencies 325 Wolberg: Image Processing Course NotesVarying GHPF Cutoff Frequencies 326 Wolberg: Image Processing Course NotesExample 327 Wolberg: Image Processing Course NotesHomomorphic Filtering (1) • Consider an image to be expressed as: f(x,y) = i(x,y) r(x,y) illumination reflectance • Take the logarithm of f(x,y) to separate i(x,y) and r(x,y): ln f(x,y) = ln i(x,y) + ln r(x,y) • Illumination component: slow spatial variations • Reflectance component: tends to vary abruptly • Filter out low frequencies to remove illumination effects. 328 Wolberg: Image Processing Course NotesHomomorphic Filtering (2) 329 Wolberg: Image Processing Course NotesExample • Interior details are obscured by glare from outside walls. • Reduction of dynamic range in brightness, together with an increase in contrast, brought out interior details and balanced the graylevels of the outside wall.  0.5 L  2.0 H 330 Wolberg: Image Processing Course NotesCorrelation • Identical to convolution except that kernel is not flipped. • Kernel is now referred to as the template. • It is the pattern that we seek to find in the larger image. f (x, y) h(x, y) F (u,v)H(u,v) 331 Wolberg: Image Processing Course NotesExample (1) 332 Wolberg: Image Processing Course NotesExample (2) Reference image Template Correlation image Peak corresponds to best match location 333 Wolberg: Image Processing Course NotesSampling Theory Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we describe sampling theory: Motivation Analysis in the spatial and frequency domains Ideal and nonideal solutions Aliasing and antialiasing Example: oversampling CD players 335 Wolberg: Image Processing Course NotesSAMPLING Sampling Theory Continuous signal g(x) Discrete (sampled) signal g (x) s Sampling theory addresses the two following problems: 1. Are the samples of g (x) sufficient to exactly describe g(x) s 2. If so, how can g(x) be reconstructed from g (x) s 336 Wolberg: Image Processing Course Notes RECONSTRUCTIONMotivation • Improper consideration leads to jagged edges in magnification and Moire effects in minification. Unfiltered image Filtered image 337 Wolberg: Image Processing Course NotesInsight: Analysis in Frequency Domain Spatial domain Frequency domain  i2fx g(x) G( f ) g(x)e dx   f f max max X   S( f ) f ( f nf )  s s s(x) (x nT ) T s s f =1/T  s s   = G ( f ) f G( f nf ) s s s  g (x)= g(x)s(x) s 2f f f 0 f f 2f s s max s max s 338 Wolberg: Image Processing Course NotesExploit WellKnown Properties of Fourier Transforms 1. Multiplication in spatial domain  convolution in frequency domain. 2. Fourier transform of an impulse train is itself an impulse train. 3. G (f) is the original spectrum G(f) replicated in the s frequency domain with period f . s G (f) s 2f f 0 f f f 2f s s max s max s G (f) G (f) high high (Base band) G (f) = G(f) + G (f)  introduced by sampling s high 339 Wolberg: Image Processing Course NotesSolution to Reconstruction Discard replicas G (f), leaving baseband G(f) intact. high Two Provisions for exact reconstruction: • The signal must be bandlimited. Otherwise, the replicas would overlap with the baseband spectrum. • f 2f (Nyquist frequency) s max 340 Wolberg: Image Processing Course NotesIdeal Filter in the Frequency Domain: Box Filter • Assuming provisions hold, we can discard G by multiplying high G (f) with a box filter in the frequency domain. s 1 f f (passband) G (f) max s H ( f )  1 0 f f (stopband) max  0 f f max max G (f) G (f) high high (Passband) 341 Wolberg: Image Processing Course NotesIdeal Filter in the Spatial Domain: Sinc Function The ideal lowpass filter in the spatial domain is the inverse Fourier transform of the ideal box: sin(x) sinc(x) x 342 Wolberg: Image Processing Course NotesReciprocal Relationship • Reciprocal relationship between spatial and frequency domains: A H(f) A/2w w w 1/2w 2/2w 3/2w • For interpolation, A must equal 1 (unity gain when sinc is centered at samples; pass through 0 at all other samples) W=f =0.5 cycle / pixel for all digital images max Highest freq: onoff sequence 1 cycle is 2 pixels (black, white, black, white), therefore ½ cycle per pixel is f . max 343 Wolberg: Image Processing Course NotesBlurring • To blur an image, W and A so that A/2W remains at 1. A 1 A 2 1/2w 1/2w 1 2 A A 2 1 W W 2 1 • To interpolate among sparser samples, A=1 and W which means that (A/2W). Since sampling decreases the amplitude of the spectrum, (A/2W)  serves to restore it. g(x) sinc(x) g (x) s   sinc()g (x)d s   Note : sinc requires infinite number of neighbours : impossible. • Ideal lowpass filtering is impossible in the spatial domain. It is not impossible in the frequency domain for finite data streams. 344 Wolberg: Image Processing Course NotesNonideal Reconstruction • Possible solution to reconstruction in the spatial domain: truncated sinc. Ripples; Gibbs phenomenon • Alternate solution: nonideal reconstruction G (f) s H(f) Passband • Nonideal because it attenuates the higher frequencies in the passband and doesn’t fully suppress the stopband. 345 Wolberg: Image Processing Course NotesAliasing If f 2f our signal is undersampled. s max If f is infinite, our signal is not bandlimited. max Either way, we have aliasing: high frequencies masquerade, or alias as, low frequencies G (f) s G (f) s adequate sampling inadequate sampling (overlapping spectrum replicas) Temporal aliasing: stagecoach wheels go backwards in films (sampled at 24 frames/sec). 346 Wolberg: Image Processing Course NotesAntialiasing Two approaches to antialiasing for combating aliasing: 1. Sample input at higher rate: increase f s G (f) G (f) s s f f f f s s s s low f (overlaps) higher f (no overlaps) s s 2. bandlimit input before sampling at low f s G (f) s G (f) s bandlimit sampling f f s s 347 Wolberg: Image Processing Course NotesAntialiasing Intuition for increasing f : s Higher sampling rates permit sloppier reconstruction filters. G (f) G (f) s s H (f) r f f f f s s s s high f allows nonideal reconstruction filter. s Marginally acceptable f s This adequate because H (f) doesn’t r Requires ideal lowpass filter taper off until FF max Pushed to limits; use sinc Many samples: linear interpolation adequate 348 Wolberg: Image Processing Course NotesNonideal Reconstruction Filter H (f) r  1 0 f f max  H ( f ) 0 f f f f  r max s max  0 f f f s max  f s f max 349 Wolberg: Image Processing Course NotesNonideal Reconstruction Filter Intuition for bandlimiting input: Inability to change fs: e.g., restricted by display resolution. For example, an image in perspective may require infinite pixels to capture all visual details. However, only finite pixels available, e.g. 512 x 512 Many pixels in oblique image map to one output pixel. Must blur neighborhood before returning a single value to the output image. In general, dull signals don’t need bandlimiting, rich ones do. 350 Wolberg: Image Processing Course NotesExample: Oversampling CD players f = 44.1 kHz (sample/second) s = 20 kHz 2 + 4.1 kHz = 44.1KHz f Nyq rate cushion max Oversampling No Oversampling Oversampling Sloppy reconstruction by Accurate digital Analog output filter Filter in CD (speaker) Player reconstructs well Large Space large overshoot/undershoots 0 f f s s Sloppy analog filter less more intermediate samples 351 Wolberg: Image Processing Course NotesAdequate Sampling Inadequate Sampling 352 Wolberg: Image Processing Course NotesAntialiasing 353 Wolberg: Image Processing Course NotesImage Resampling Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review image resampling: Ideal resampling Mathematical formulation Resampling filter Tradeoffs between accuracy and complexity Software implementation 355 Wolberg: Image Processing Course NotesDefinition • Image Resampling: Process of transformation a sampled image from one coordinate system to another. 356 Wolberg: Image Processing Course NotesIdeal Resampling • Ideal Resampling: reconstruction, warping, prefiltering, sampling 357 Wolberg: Image Processing Course NotesMathematical Formulation (1) Stages : Math Definition : Discrete Input f (u),u z Reconstruction Input f f (u) r(u) f (k)r(u k) c kz 1 Warped Signal g (x) f (m (x)) c c ' Continuous Output g (x) g (x) h(x) g (t)h(x t)dt c c c  ' Discrete Output g(x) g (x)S(x) c Two filtering components : reconstruction and prefiltering (bandlimiting warped signal before sampling). Cascade them into one by working backwards from g(x) to f (u) : ' g(x) g (x) for x z c  11 g(x) f (m (t))h(x t)dt f (k)r(m (t) k) h(x t)dt  c   kz 1 g(x) f (k)(x, k) where (x, k) r(m (t) k)h(x t)dt   kz 358 Wolberg: Image Processing Course NotesMathematical Formulation (2) 1 Spatially varying resampling (x,k) r(m (t) k)h(x t)dt  filter expressed in terms of output space Selects filter response used to index filter coefficients (Wide) (Narrow) 359 Wolberg: Image Processing Course NotesMathematical Formulation (3) • We can express (x,k) in terms of input space: Let t= m(u), we have m (x, k) r(u k)h(x m(u)) dt  u m where is the determinant of Jacobian matrix : u x x m dmmx u v 1D : 2D : where x u y y u duuu u v 360 Wolberg: Image Processing Course NotesResampling Filter Linear warps (space invariant) : '1  (m (x k)) ' (x, k) h (u) r(u) J h(uJ) r(u)  (x,k) = r(u): Shape of r(u) remains the same, independently of m(u), mag (independently of scale factor)  (x,k) = Jh(uJ): Shape of prefilter is based on desired frequency response min characteristic (performance in passband and stopband). Unlike r(u), though, the prefilter must be scaled proportional to the minification factor (broader and shorter for more minification) 361 Wolberg: Image Processing Course NotesFourier Transform Pairs • The shape of  is a direct consequence of the reciprocal relation between min the spatial and frequency domains.  denotes Fourier Transform pair i2fu H (u) h(u)e du  i2fu H (m(u)) h(m(u))e du  m Let x au m(u) and dx du u 1 dx xm i2fm ( x)1 H (m(u)) h(x)e where m (x) ;  J  a  m au u fx i2 1 a h(au) h(x)e dx  a 1 f  h(au) H  a a  362 Wolberg: Image Processing Course NotesReciprocal Relationship H(f) h(x) h(ax) , a1 Magnification yields less highfreq (more blurred) h(ax) , a1 Minification introduces higher freqs (more visual detail) Intuition: f = 1 / T Consequences: narrow filters in spatial domain (desirable) yield wide frequency spectrums (undesirable). Tradeoff between accuracy and complexity. 363 Wolberg: Image Processing Course NotesSoftware (1) resample1D(IN, OUT, INlen, OUTlen, filtertype, offset) unsigned char IN, OUT; int INlen, OUTlen, filtertype, offset; int i; int left, right; // kernel extent in input int pixel; // input pixel value double u, x; // input (u) , output (x) double scale; // resampling scale factor double (filter)(); // pointer to filter fct double fwidth; // filter width (support) double fscale; // filter amplitude scale double weight; // kernel weight double acc; // convolution accumulator scale = (double) OUTlen / INlen; 364 Wolberg: Image Processing Course NotesSoftware (2) switch(filtertype) case 0: filter = boxFilter; // box filter (nearest nbr) fwidth = .5; break; case 1: filter = triFilter; // triangle filter (lin intrp) fwidth = 1; break; case 2: filter = cubicConv; // cubic convolution filter fwidth = 2; break; case 3: filter = lanczos3; // Lanczos3 windowed sinc fct fwidth = 3; break; case 4: filter = hann4; // Hann windowed sinc function fwidth = 4; // 8point kernel break; 365 Wolberg: Image Processing Course NotesSoftware (3) if(scale 1.0) // minification: h(x) h(xscale)scale fwidth = fwidth / scale; // broaden filter fscale = scale; // lower amplitude / roundoff fwidth to int to avoid intensity modulation / if(filtertype == 0) fwidth = CEILING(fwidth); fscale = 1.0 / (2fwidth); else fscale = 1.0; // project each output pixel to input, center kernel,and convolve for(x=0; xOUTlen; x++) / map output x to input u: inverse mapping / u = x / scale; / left and right extent of kernel centered at u / if(u fwidth 0) left = FLOOR (u fwidth); else left = CEILING(u fwidth); right = FLOOR(u + fwidth); 366 Wolberg: Image Processing Course NotesSoftware (4) / reset acc for collecting convolution products / acc = 0; / weigh input pixels around u with kernel / for(i=left; i = right; i++) pixel = IN CLAMP(i, 0, INlen1)offset; weight = (filter)((u i) fscale); acc += (pixel weight); / assign weighted accumulator to OUT / OUTxoffset = acc fscale; 367 Wolberg: Image Processing Course NotesSoftware (5) double boxFilter(double t) if((t .5) (t = .5)) return(1.0); return(0.0); double triFilter(double t) if(t 0) t = t; if(t 1.0) return(1.0 t); return(0.0); 368 Wolberg: Image Processing Course NotesSoftware (6) double cubicConv(double t) double A, t2, t3; if(t 0) t = t; t2 = t t; t3 = t2 t; A = 1.0; // userspecified free parameter if(t 1.0) return((A+2)t3 (A+3)t2 + 1); if(t 2.0) return(A(t3 5t2 + 8t 4)); return(0.0); 369 Wolberg: Image Processing Course NotesSoftware (7) double sinc(double t) t = PI; if(t = 0) return(sin(t) / t); return(1.0); double lanczos3(double t) if(t 0) t = t; if(t 3.0) return(sinc(t) sinc(t/3.0)); return(0.0); 370 Wolberg: Image Processing Course NotesSoftware (8) double hann4(double t) int N = 4; // fixed filter width if(t 0) t = t; if(t N) return(sinc(t) (.5 + .5cos(PIt / N))); return(0.0); 371 Wolberg: Image Processing Course NotesImage Reconstruction Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we describe image reconstruction: Interpolation as convolution Interpolation kernels for: • Nearest neighbor • Triangle filter • Cubic convolution • BSpline interpolation • Windowed sinc functions 373 Wolberg: Image Processing Course NotesIntroduction • Reconstruction is synonymous with interpolation. • Determine value at position lying between samples. • Strategy: fit a continuous function through the discrete input samples and evaluate at any desired set of points. • Sampling generates infinite bandwidth signal. • Interpolation reconstructs signal by smoothing samples with an interpolation function (kernel). 374 Wolberg: Image Processing Course NotesInterpolation For equispaced data, interpolation can be expressed as a convolution: K1 f (x) c h(x x ) where K is the number of neighborhood pixels  k r k0 and c are coefficients for kernel h k Kernel Samples: h(d), h(1d), h(1d), h(2d) If h is symmetric, h(d), h(1+d), h(1d), h(2d) 375 Wolberg: Image Processing Course NotesInterpolation Kernel • Set of weights applied to neighborhood pixels • Often defined analytically • Usually symmetric: h(x) = h(x) • Commonly used kernels: Nearest neighbor (pixel replication) Triangle filter (linear interpolation) Cubic convolution (smooth; used in digital cameras) Windowed sinc functions (highest quality, more costly) 376 Wolberg: Image Processing Course NotesNearest Neighbor x x x x k1 k k k1 Interpolating Polynomial : f (x) f (x )  x k 2 2 1 0 x 0.5 Interpolating Kernel : h(x)  0 0.5 x  x x x k1 k k+1 Other names: box filter, sampleand hold function, and Fourier window. Poor stopband. NN achieves magnification by pixel replication. Very blocky. Shift errors of up to 1/2 pixel are possible. Common in hardware zooms. 377 Wolberg: Image Processing Course NotesTriangle Filter Interpolating Polynomial : f (x) a x a 1 0 x x  0 1 f f  a a Solve for a , a 0 1 1 0 1 0  1 1  x x x k1 k k+1  x x 0 f (x) f ( f f ) 0 1 0  x x 1 0  1 x 0 x1 Interpolation Kernel : h(x)  0 1 x  Other names for h: triangle filter, tent filter, roof function, chateau function, and Bartlett window. In the frequency domain: = 2 Sinc x Sinc = Sinc 378 Wolberg: Image Processing Course NotesCubic Convolution (1) Third degree approximation to sinc. Its kernel is derived from constraints imposed on the general cubic spline interpolation formula. 3 2  a x a x a x a 0 x 1 30 20 10 00   3 2 h(x) a x a x a x a 1 x 2  31 21 11 01  0 2 x   Determine coefficients by applying following constraints: 1. h(0) = 1 and h(x) = 0 for x = 1, 2 2. h must be continuous at x = 0, 1, 2 3. h must have a continuous first derivative at x = 0, 1, 2 379 Wolberg: Image Processing Course NotesCubic Convolution (2) Constraint (1) states that when h is centered on an input sample, the interpolation function is independent of neighboring samples. First 2 constraints give 4 equations: 1 h(0) a 00  0 h(1 ) a a a a 30 20 10 00  0 h(1 ) a a a a 31 21 11 01  0 h(2 ) 8a 4a 2a a 31 21 11 01 3 more equations are obtained from constraint (3) : ' '  a h (0 ) h (0 ) a 10 10 ' ' 3a 2a a h (1 ) h (1 ) 3a 2a a 30 20 10 31 21 11 ' ' 12a 4a a h (2 ) h (2 ) 0 31 21 11 Total :7 equations, 8 unknowns  free variable (a a ) 31 3 2  (a 2) x (a 3) x1 0 x1   3 2 h(x) a x 5a x 8a x 4a 1 x 2   0 2 x   380 Wolberg: Image Processing Course NotesCubic Convolution (3) How to pick a Add some heuristics (make it resemble Sinc function): ’’ h (0) = 2(a+3) 0  a 3 Concave downward at x = 0 ’’ h (1) = 4a 0 Concave upward at x = 1 This bounds a to the 3, 0 range. Common choices: a = 1 matches the slope of sinc at x=1 (sharpens image) a = 0.5 makes the Taylor series approximately agree in as many terms as possible with the original signal a = .75 sets the second derivative of the 2 cubic polynomials in h to 1 nd (continuous 2 derivative at x = 1) 381 Wolberg: Image Processing Course NotesCubic Splines (1) 3 2 f (x) a (x x ) a (x x ) a (x x ) a k 3 k 2 k 1 k 0 rd 6 polynomial segments, each of 3 degree. ’ ’ “ f s are joined at x (for k=1,…, n2) such that f , f , and f are continuous. k k k k k a y 0 k ' a y where y y y 1 k k1 k k ' ' a 3y 2y y 2 k k k1 ' ' a2y y y 3 k k k1 (proof in App. 2) 382 Wolberg: Image Processing Course NotesCubic Splines (2) • The derivatives may be determined by solving the following system of linear equations: ' 2 4 5y 4y y y  0 1 2 0    ' 1 4 1 3(y y ) y 2 0 1    '    1 4 1 3(y y ) y 3 1 2    ' 1 4 1 3(y y ) y 4 2  3         ' 1 4 1 3(y y 3) y   n1 n n2  '  4 2 y 4y 5y y  n3 n2 n1 n1   global dependencies • Introduced by the constraints for continuity in the first and second derivatives at the knots. 383 Wolberg: Image Processing Course NotesBSplines To analyze cubic splines, introduce cubic BSpline interpolation kernel: B box filter 0 B = B B B 2 0 0 0 B = B B 1 0 0 B = B B B B 3 0 0 0 0 3 2  3 x 6 x 4 0 x1  1  3 2 Parzen Window: Not interpolatory because it h(x) x 6 x12 x 8 1 x 2  does not satisfy h(0) = 1, and h(1) = h(2) = 0. 6  0 2 x  Indeed, it approximates the data.  384 Wolberg: Image Processing Course NotesInterpolatory BSplines j2 f (x ) c h(x x ) j k j k k j2 4 1 Since h(0) , h(1) h(1) , h(2) h(2) 0 6 6 1 we have f (x ) (c 4c c ) j j1 j j1 6 f 4 1 c  0 0  f 1 4 1 c 1 1   f 1 4 1 c 1 2 2     6   f 1 4 1 c n2 n2  f 1 4 c    n1 n1 F K C 1 C K F 1 K inverse of tridiagonal matrix; Computation is O(n) 1 All previous methods used data values for c from C = K F. k 385 Wolberg: Image Processing Course NotesTruncated Sinc Function • Alternative to previous kernels: use windowed sinc function. 1 Truncated Sinc X 1 0 x 0.5 Rect(x)  0 0.5 x  Truncating in spatial domain = convolving spectrum (box) with a Sinc function. Ringing can be mitigated by using a smoothly tapering windowing function. Popular window functions: Hann, Hamming, Blackman, Kaiser, and Lanczos. 386 Wolberg: Image Processing Course NotesHann/Hamming Window 2x N1   (1)cos x Hann/Hamming(x)  N1 2  0 o / w  N = number of samples in windowing function. Hann:  = 0.5; Hamming:  = 0.54. Also known as raised cosine window. h(x) (sinc(x))(Hann(x)) H(f) is sinc+2 shifted sincs. These cancel the right and left side lobes of Rect(x). 387 Wolberg: Image Processing Course NotesBlackman Window 2x 4x N1   .42.5cos.08cos x Blackman(x)  N1 N1 2  0 o / w  388 Wolberg: Image Processing Course NotesLanczos Window (1) h(x) = sinc(x) Lanczos2(x) x  sinc 0 x 2  Lanczos2(x) sinc(x) sinc(x/2) Rect(x/4) Spatial Domain 2    0 2 x  window = x 10 8 6 4 2 0 2 4 6 8 10 10 8 6 4 2 0 2 4 6 8 10 Rect(x/4) sinc(x/2) h(x)=sinc(x)Lanczos2(x) H(f) = Rect(f) Rect(2f) sinc(4f) Frequency Domain 389 Wolberg: Image Processing Course NotesLanczos Window (2) • Generalization to N lobes: x  sinc 0 x N  LanczosN(x) N    0 N x  • Let N = 3, this lets 3 lobes pass under the Lanczos window. • Better passband and stopband response 390 Wolberg: Image Processing Course NotesComparison of Interpolation Methods NN, linear, cubic convolution, windowed sinc, sinc poor …………………………………………….. ideal (blocky, blurred, ringing, no artifacts) 391 Wolberg: Image Processing Course NotesConvolution Implementation 1. Position (center) kernel in input. 2. Evaluate kernel values at positions coinciding with neighbors. 3. Compute products of kernel values and neighbors. 4. Add products; init output pixel. Step (1) can be simplified by incremental computation for spaceinvariant warps. (newpos = oldpos + inc). Step (2) can be simplified by LUT. 392 Wolberg: Image Processing Course NotesInterpolation with Coefficient Bins Implementation 1: Interp. with Coefficient Bins (for spaceinvariant warps) • Strategy: accelerate resampling by precomputing the input weights and storing them in LUT for fast access during convolution. n intervals (bins) between input pixels Input x Output u  bin: (quantize u) h , h h , h Kernel values 2 1 1 2 ……………… n1 i 0 1 2 Interleaved f f 2 1 f Let d = 1 – i/n (ud1) 1 f h 2 1 h = h(d); h = h(1d) h 1 1 1 h 2 h h = h(1+d); h = h(2d) 2 2 2 393 Wolberg: Image Processing Course NotesUninterleaved Coefficient Bins ……………… ……………… ……………… n1 n1 n1 i i i 0 1 2 0 1 2 0 1 2 Oversample Oversample Oversample ii = bin 0  ii oversample val = IN2 kern2 oversample + ii + IN1 kern1 oversample + ii + IN 0 kernii Refer to code on p. 151 + IN 1 kern1 oversample ii + IN 2 kern2 oversample ii + IN 3 kern3 oversample ii; if(ii == 0)val +=IN3 kern3 oversample ii; kern2+ii Kern1ii Kernii kern2+ii Kern2ii Kern3ii 0  ii  1 kern3+ii IN3 IN2 IN1 IN0 IN1 IN2 IN3 394 Wolberg: Image Processing Course NotesKernel Position • Since we are assuming space invariance, the new position for the kernel = oldpos + offset. INlenoversample  offset dii dff ; dii whole bins  OUTlen  dff partial bin (INlenoversample)OUTlen Subpixel (bin) Offset must be accurate to avoid accrual of error in the incremental repositioning Pixel Pixel of the kernel. Blowup J K L I of bins. Current pos. Let dii = 1 dff = .6 I J K L Old New pos. pos. Offset = 1.6 bins 395 Wolberg: Image Processing Course NotesForward vs. Inverse Mapping v Y Input Output u X Ch. 3, Sec. 1 • Forward mapping: x = X(u, v); y = Y(u, v) • Inverse mapping: u = U(x,y); v = V(x, y) Output (accumulator) Input Input Output Inverse mapping Forward mapping Coefficient Bins for kernel eval for fast Convolution for image reconstruction. 0 2OV OV 3OV 1d d LHS RHS OV OV dOV OV 396 Wolberg: Image Processing Course NotesFant’s Algorithm Implementation 2: Fant’s Resampling Algorithm (for spacevar. warps) Resampler Input and output are streams of pixels that are consumed and generated at rate determined by the spatial mapping. Three conditions per stream: 1) Current input pixel is entirely consumed without completing an output pixel. 2) The input is entirely consumed while completing the output pixel. 3) Output pixel computed without entirely consuming the current input pixel. Algorithm uses linear interpolation for image reconstruction and box filtering (unweighted averaging) for antialiasing. Code on p.156. 397 Wolberg: Image Processing Course NotesExample A (100)(.4) 40 0 m(u) 2.3 3.9 0.6 3.2 3.3  .4 .4  A (100) 1 (106) (1) 101  1  1.7 1.7  Input 100 106 92 90  1.4 1.4  A (100) 1 (106) (.3) (106)(.7) 106  2  .6 2.3 3.2 3.3 3.9 1.7 1.7   Output  .7 .7  A (106) 1 (92) (.2) (92)(.1) (90)(.6) 82  0 2 3 4 3 1 A A A A .9 .9 0 1 2 3  106 100 92 90 398 Wolberg: Image Processing Course NotesAntialiasing Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review antialiasing: Supersampling (uniform, adaptive, irregular) Direct convolution • FeibushLevoyCook • GangnetPerryCoveignoux • GreeneHeckbert (Elliptical Weighted Average) Prefiltering • Mipmaps / pyramids • Summed area tables 400 Wolberg: Image Processing Course NotesAntialiasing • Aliasing is due to undersampling the input signal. • It results in false artifacts, e.g., moire effects. • Antialiasing combats aliasing. Solutions: • Raise sampling rate • Bandlimit input • In practice, do both in inverse mapping formulation • Implementations: • Supersampling (uniform, adaptive, irregular) • Direct convolution • Prefiltering 401 Wolberg: Image Processing Course NotesPoint Sampling • One input sample per output pixel • Problem: information is lost between samples • Solution: sample more densely and average results 1D: 2D: 402 Wolberg: Image Processing Course NotesSupersampling (1) • Multiple input samples per output pixel 1D: 2D: LPF 403 Wolberg: Image Processing Course NotesSupersampling (2) 4 samples/pixel 1 sample/pixel 64 samples/pixel 256 samples/pixel 16 samples/pixel 404 Wolberg: Image Processing Course NotesAdaptive Supersampling • Collect more samples in areas of high intensity variance or contrast. A B If (f(A,B,C,D,E)thr) then subdivide / quad tree / E C D • To avoid artifacts due to regular sampling pattern (along rectilinear grid), use irregular sampling. • Three common forms of stochastic sampling: Poisson Jittered Pointdiffusion sampling 405 Wolberg: Image Processing Course NotesDirect Convolution: FeibushLevoyCook (1980) Inverse Mapping Input bounding box of kernel in output (screen) space 1. All input pixels contained within the bounding rect. of this quadrilateral are mapped onto output. 2. The extra selected points are eliminated by clipping them against the bounding rect. of kernel. 3. Init output pixel with weighted average of the remaining samples. Note that filter weights are stored in a LUT (e.g., a Gaussian filter.) Advantages: 1) Easy to collect input pixels in rect. input region 2) Easy to clip output pixels in rect. output region. 3) Arbitrary filter weights can be stored in LUT. 406 Wolberg: Image Processing Course NotesDirect Convolution: GangnetPerryCoveignoux (1982) Improve results by supersampling. Major axis determines supersampling rate. They used bilinear interpolation for image reconstruction and a truncated sinc (2 pixels on each side) as kernel. 407 Wolberg: Image Processing Course NotesElliptical weighted average (EWA): GreeneHeckbert (1986) Q(u,v) EWA distorts the circular kernel into an ellipse in the input where the weighting can be computed directly. No mapping back to output space. 2 2 Q(u,v) Au Buv Cv u v 0 is the center of the ellipse 2 2 A VV x y B2(U VU V ) x x y y 2 2 C UU x y where  uvuv   (U ,V ) , (U ,V ) ,  x x y y  xxyy   2 Point inclusion test : Q(u,v) F for F (U VU V ) x y y x 408 Wolberg: Image Processing Course NotesEWA Advantages • Only 1 function evaluated; not 4 needed for quadrilateral. • If QF, then sample is weighted with appropriate LUT entry. • In output space, LUT is indexed by r • Instead of determining which concentric circle the point lies in the output, we determine which concentric ellipse the point lies in the input and use it to index the LUT. r Q r 409 Wolberg: Image Processing Course NotesComparisons • All direct convolution methods are O(N) where N = input pixels accessed. In Feibush and Gangnet, these samples must be mapped into output space, EWA avoids this costly step. • Direct convolution imposes minimal constraints on filter area (quad, ellipse) and filter kernel (precomputed LUT). Additional speedups: prefiltering with pyramids and preintegrated tables approx convolution integral w/ a constant of accesses (indep. of input pixels in preimage.) • Drawback of pyramids and preintegrated tables: filter area must be square or rectangular; kernel must be a box filter. 410 Wolberg: Image Processing Course NotesPyramids d ¼ res. v ½ res. 1 1 1 4 memory cost  1 4 16 64 3 Original only 33 more memory (full res.) u Size of square is used to determine pyramid level to sample. Aliasing vs. Blurring tradeoff. 2 2 2 2  uvuv  2  d max ,   xxyy    d is proportional to span of the preimage area. 411 Wolberg: Image Processing Course NotesSummedArea Table y 1 Table stores running sum R R 2 y 0 sum Tx , y Tx , y Tx , y  Tx , y R 1 1 1 0 0 1 0 0 R R 0 1 Restricted to rectangular regions and box filtering x x 0 1 To compute current entry Tx1 , y1 let sum R be current pixel value v: Tx , y  v Tx , y  Tx , y Tx , y 1 1 1 0 0 1 0 0 412 Wolberg: Image Processing Course NotesExample 90 10 20 30 265 460 820 1190 175 360 700 1040 50 60 70 80 125 250 520 780 25 75 200 180 100 150 220 300 100 50 70 80 Input Summedarea table sum Tx , y Tx , y Tx , y  Tx , y R 1 1 1 0 0 1 0 0 Tx , y  v Tx , y  Tx , y Tx , y 1 1 1 0 0 1 0 0 413 Wolberg: Image Processing Course NotesSpatial Transformations Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review spatial transformations: Forward and inverse mappings Transformations • Linear • Affine • Perspective • Bilinear Inferring affine and perspective transformations 415 Wolberg: Image Processing Course NotesForward and Inverse Mappings • A spatial transformation defines a geometric relationship between each point in the input and output images. • Forward mapping: x, y = X(u,v), Y(u,v) • Inverse mapping : u, v = U(x,y), V(x,y) • Forward mapping specifies output coordinates (x,y) for each input point (u,v). • Inverse mapping specifies input point (u,v) for each output point (x,y). Output (accumulator) Input Input Output Forward mapping Inverse mapping 416 Wolberg: Image Processing Course NotesLinear Transformations a a  11 12 x, y u,v  a a  21 22 x a u a v 11 21 y a u a v 12 22 • Above equations are linear because they satisfy the following two conditions necessary for any linear function L(x): 1) L(x+y) = L(x) + L(y) 2) L(cx) = cL(x) for any scalar c and position vectors x and y. • Note that linear transformation are a sum of scaled input coordinate: they do not account for simple translation. We want: a a  11 12 x a u a v a  11 21 31 x, y u,v,1 a a 21 22  y a u a v a 12 22 32  a a  31 32 417 Wolberg: Image Processing Course NotesHomogeneous Coordinates • To compute inverse, the transformation matrix must be square. • Therefore, a a 0  11 12  x, y,1 u,v,1 a a 0 21 22  a a 1  31 32 • All 2D position vectors are now represented with three components: homogeneous notation. • Third component (w) refers to the plane upon which the transformation operates. • u, v, 1 = u, v position vector lying on w=1 plane. • The representation of a point in the homogeneous notation is no longer unique: 8, 16, 2 = 4, 8, 1 = 16, 32, 4. ’ ’ ’ • To recover any 2D position vector px,y from p =x, y, w, divide by the h ’. homogeneous coordinate w ' '  x y x, y  ' ' w w  418 Wolberg: Image Processing Course NotesAffine Transformations (1) a a 0  11 12  x, y,1 u,v,1 a a 0 21 22   a a 1  31 32 1 0 0   Translation  0 1 0 1 0 0     T T 1 Shear rows  H 1 0 u v  u   0 0 1  cos sin 0   1 H 0  u Rotation  sin cos 0   Shear columns  0 1 0  0 0 1    0 0 1  S 0 0  u  Scale  0 S 0 v   0 0 1  419 Wolberg: Image Processing Course NotesAffine Transformations (2) • Affine transformation have 6 degrees of freedom: a , a , a , a , a , a . 11 21 31 12 22 23 • They can be inferred by giving the correspondence of three 2D points between the input and output images. That is, (u ,v ) (x , y )  1 1 1 1  (u ,v ) (x , y ) 6 constraints : ( 3 for u x, 3 for v y)  2 2 2 2  (u ,v ) (x , y ) 3 3 3 3  • All points lie on the same plane. • Affine transformations map triangles onto triangles. 420 Wolberg: Image Processing Course NotesAffine Transformations (3) Skew (shear) Rotation/Scale Rotation 421 Wolberg: Image Processing Course NotesPerspective Transformations (1) ' x a u a v a a a a 11 21 31  11 12 13 x ' ' ' ' w a u a v a 13 23 33 x , y , w  u,v, w a a a 21 22 23  '  y a u a v a a a a 12 22 32  31 32 33 y ' w a u a v a 13 23 33 a 0  13  a not necessarily 0 as in affine transformations 23   a 1  33 • Without loss of generality, we set a =1. 33 • This yields 8 degrees of freedom and allows us to map planar quadrilaterals to planar quadrilaterals (correspondence among 4 sets of 2D points yields 8 coordinates). • Perspective transformations introduces foreshortening effects. • Straight lines are preserved. 422 Wolberg: Image Processing Course NotesPerspective Transformations (2) 423 Wolberg: Image Processing Course NotesBilinear Transforms (1) a b  3 3  a b 2 2  x, y uv,u,v,1  a b 1 1  a b  0 0 nd • 4corner mapping among nonplanar quadrilaterals (2 degree due to uv factors). • Conveniently computed using separability (see p. 5760). • Preserves spacing along edges. • Straight lines in the interior no longer remain straight. 424 Wolberg: Image Processing Course NotesBilinear Transforms (2) 425 Wolberg: Image Processing Course NotesExamples Similarity transformation (RST) Affine transformation Perspective transformation Polynomial transformation 426 Wolberg: Image Processing Course NotesScanline Algorithms Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review scanline algorithms: Incremental texture mapping 2pass CatmullSmith algorithm • Rotation • Perspective 3pass shear transformation for rotation Morphing 428 Wolberg: Image Processing Course NotesCatmullSmith Algorithm • Twopass transform • First pass resamples all rows: u, v  x, v x, v = F (u), v where F (u) = X(u, v) is the forward mapping fct v v • Second pass resamples all columns: x, v  x, y x, y = x, G (v) where G (v)=Y(H (v), v) x x x • H (v) is the inverse projection of x’, the column we wish x to resample. • It brings us back from x, v to u, v so that we can directly index into Y to get the destination y coordinates. 429 Wolberg: Image Processing Course NotesExample: Rotation (1) cos sin  x, y u,v   sin cos  Pass 1: x,v u cos vsin,v Pass 2 : a) Compute H (v). Recall that x u cos vsin x x vsin u cos b) Compute G (v). Substitute H (v) into y u sin v cos x x xsin v y cos 430 Wolberg: Image Processing Course Notes2Pass Rotation scale/shear rows scale/shear columns 431 Wolberg: Image Processing Course Notes3Pass Rotation Algorithm • Rotation can be decomposed into two scale/shear matrices. 1 tan  cos sin cos 0  1  R  0  sin cos sin 1    cos • Three pass transform uses on shear matrices. 1 0 1 0  cos sin 1 sin    R   tan 1 tan 1   sin cos 0 1   2 2   • Advantage of 3pass: no scaling necessary in any pass. 432 Wolberg: Image Processing Course Notes3Pass Rotation shear rows shear columns shear rows 433 Wolberg: Image Processing Course NotesSoftware Implementation • The following slides contain code for initMatrix.c to produce a 3x3 perspective transformation matrix from a list of four corresponding points (e.g. image corners). • That matrix is then used in perspective.c to resample the image. • The code in resample.c performs the actual scanline resample. 434 Wolberg: Image Processing Course NotesinitMatrix.c (1) / initMatrix: Given Icorr, a list of the 4 correspondence points for the corners of image I, compute the 3x3 perspective matrix in Imatrix. / void initMatrix(imageP I, imageP Icorr, imageP Imatrix) int w, h; float p, a, a13, a23; float x0, x1, x2, x3; float y0, y1, y2, y3; float dx1, dx2, dx3, dy1, dy2, dy3; / init pointers / a = (float ) Imatrixbuf0; p = (float ) Icorrbuf0; / init u,v,x,y vars and print them / x0 = p++; y0 = p++; x1 = p++; y1 = p++; x2 = p++; y2 = p++; x3 = p++; y3 = p++; 435 Wolberg: Image Processing Course NotesinitMatrix.c (2) w = Iwidth; h = Iheight; UIprintf("\nCorrespondence points:\n"); UIprintf("4d 4d 6.1f 6.1f\n", 0, 0, x0, y0); UIprintf("4d 4d 6.1f 6.1f\n", w, 0, x1, y1); UIprintf("4d 4d 6.1f 6.1f\n", w, h, x2, y2); UIprintf("4d 4d 6.1f 6.1f\n", 0, h, x3, y3); / compute auxiliary vars / dx1 = x1 x2; dx2 = x3 x2; dx3 = x0 x1 + x2 x3; dy1 = y1 y2; dy2 = y3 y2; dy3 = y0 y1 + y2 y3; 436 Wolberg: Image Processing Course NotesinitMatrix.c (3) / compute 3x3 transformation matrix: a0 a1 a2 a3 a4 a5 a6 a7 a8 / a13 = (dx3dy2 dx2dy3) / (dx1dy2 dx2dy1); a23 = (dx1dy3 dx3dy1) / (dx1dy2 dx2dy1); a0 = (x1x0+a13x1) / w; a1 = (y1y0+a13y1) / w; a2 = a13 / w; a3 = (x3x0+a23x3) / h; a4 = (y3y0+a23y3) / h; a5 = a23 / h; a6 = x0; a7 = y0; a8 = 1; 437 Wolberg: Image Processing Course Notesperspective.c (1) define X(A, U, V) ((A0U + A3V + A6) / (A2U + A5V + A8)) define Y(A, U, V) ((A1U + A4V + A7) / (A2U + A5V + A8)) define H(A, X, V) (((A5V+A8)X+ A3V + A6) / (A2X A0)) / perspective: Apply a perspective image transformation on input I1. The 3x3 perspective matrix is given in Imatrix. The output is stored in I2. / void perspective(imagP I1, imageP Imatrix, imageP I2) int i, w, h, ww, hh; uchar p1, p2; float u, v, x, y, xmin, xmax, ymin, ymax, a, F; imageP II; w = I1width; h = I1height; a = (float ) Imatrixbuf0; 438 Wolberg: Image Processing Course Notesperspective.c (2) xmin = xmax = X(a, 0, 0); x = X(a, w, 0);xmin = MIN(xmin, x); xmax = MAX(xmax, x); x = X(a, w, h);xmin = MIN(xmin, x); xmax = MAX(xmax, x); x = X(a, 0, h);xmin = MIN(xmin, x); xmax = MAX(xmax, x); ymin = ymax = Y(a, 0, 0); y = Y(a, w, 0);ymin = MIN(ymin, y); ymax = MAX(ymax, y); y = Y(a, w, h);ymin = MIN(ymin, y); ymax = MAX(ymax, y); y = Y(a, 0, h);ymin = MIN(ymin, y); ymax = MAX(ymax, y); ww = CEILING(xmax) FLOOR(xmin); hh = CEILING(ymax) FLOOR(ymin); / allocate mapping fct buffer / x = MAX(MAX(w, h), MAX(ww, hh)); F = (float ) malloc(x sizeof(float)); if(F == NULL) IPbailout("perspective: No memory"); / allocate intermediate image / II = IPallocImage(ww, h, BWTYPE); IPclearImage(II); p1 = (uchar ) I1buf0; p2 = (uchar ) IIbuf0; 439 Wolberg: Image Processing Course Notesperspective.c (3) / first pass: resample rows / for(v=0; vh; v++) / init forward mapping function F; map xmin to 0 / for(u=0; u w; u++) F(int) u = X(a, u, v) xmin; resample(p1, w, 1, F, p2); p1 += w; p2 += ww; / display intermediate image / IPcopyImage(II, NextImageP); IPdisplayImage(); / init final image / IPcopyImageHeader(I1, I2); I2width = ww; I2height = hh; IPinitChannels(I2, BWTYPE); IPclearImage(I2); 440 Wolberg: Image Processing Course Notesperspective.c (4) / second pass: resample columns / for(x=0; xww; x++) p1 = (uchar ) IIbuf0 + (int) x; p2 = (uchar ) I2buf0 + (int) x; / skip past padding / for(v=0; vh; v++,p1+=ww) if(p1) break; / check for nonzero pixel / u = H(a, (x+xmin), v); / else, if pixel is black / if(u=0 uw) break; / then check for valid u / / init forward mapping function F; map ymin to 0 / for(i=0; vh; v++) u = H(a, (x+xmin), v); u = CLIP(u, 0, w1); Fi++ = Y(a, u, v) ymin; resample(p1, i, ww, F, p2); IPfreeImage(II); 441 Wolberg: Image Processing Course Notesresample.c (1) / Resample the len elements of src (with stride offst) into dst according to the monotonic spatial mapping given in F (len entries). The scanline is assumed to have been cleared earlier. / void resample(uchar src, int len, int offst, float F, uchar dst) int u, uu, x, xx, ix0, ix1, I0, I1, pos; double x0, x1, dI; if(F0 Flen1) pos = 1; / positive output stride / else pos = 0; / negative output stride / for(u=0; ulen1; u++) / index into src / uu = u offst; / output interval (real and int) for input pixel u / if(pos) / positive stride / ix0 = x0 = Fu; ix1 = x1 = Fu+1; I0 = srcuu; I1 = srcuu+offst; 442 Wolberg: Image Processing Course Notesresample.c (2) else / flip interval to enforce positive stride / ix0 = x0 = Fu+1; ix1 = x1 = Fu; I0 = srcuu+offst; I1 = srcuu; / index into dst / xx = ix0 offst; / check if interval is embedded in one output pixel / if(ix0 == ix1) dstxx += I0 (x1x0); continue; / else, input straddles more than one output pixel / / left straddle / dstxx += I0 (ix0+1x0); 443 Wolberg: Image Processing Course Notesresample.c (3) / central interval / xx += offst; dI = (I1I0) / (x1x0); for(x=ix0+1; xix1; x++,xx+=offst) dstxx = I0 + dI(xx0); / right straddle / if(x1 = ix1) dstxx += (I0 + dI(ix1x0)) (x1ix1); 444 Wolberg: Image Processing Course NotesImage Warping: A Review Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review digital image warping: Geometric transformations Forward inverse mapping Sampling Image reconstruction Interpolation kernels Separable transforms Fant’s resampling algorithm 446 Wolberg: Image Processing Course NotesDefinition • Image warping deals with the geometric transformation of digital images. 447 Wolberg: Image Processing Course NotesGeometric Transformations • Affine • Perspective • Bilinear • Polynomial • Splines • Elastic (local deformations) 448 Wolberg: Image Processing Course NotesSpatial Transformations • Forward Mapping x, y = X(u, v), Y(u, v) • Inverse Mapping u, v = U(x, y), V(x, y) 449 Wolberg: Image Processing Course NotesForward / Inverse Mapping Input image Output (accumulator) image v y u x Input image Output image v y u x 450 Wolberg: Image Processing Course NotesSampling • Point Sampling • Area Sampling LPF 451 Wolberg: Image Processing Course NotesArea Sampling • Treats pixels as finite areas • Avoids aliasing (undersampling) artifacts • Approximated by supersampling 452 Wolberg: Image Processing Course NotesSupersampling • Average of projected subpixels LPF 1 2 4 8 16 453 Wolberg: Image Processing Course NotesImage Reconstruction • Pixel values are known at integer positions • Samples can project to realvalued positions • How do we evaluate the image values at these real valued positions Reconstruction Reconstruct 454 Wolberg: Image Processing Course NotesInterpolation • Reconstruction interpolates the input • In practice, interpolation is performed at points of interest only, not entire function • Interpolation is achieved by convolution 455 Wolberg: Image Processing Course NotesConvolution N g(x) f (x )h(x x )  k k k0 g(x) 456 Wolberg: Image Processing Course NotesInterpolation Functions Interpolation functions/kernels include: • Box filter • Triangle filter • Cubic convolution • Windowed sinc functions 457 Wolberg: Image Processing Course NotesBox Filter • Nearest neighbor interpolation • Blocky artifacts may occur 458 Wolberg: Image Processing Course NotesTriangle Filter • Linear interpolation • Popular for use with small deformations 459 Wolberg: Image Processing Course NotesCubic Convolution • Local cubic interpolation algorithm • Advanced feature in digital cameras 460 Wolberg: Image Processing Course NotesWindowed Sinc Function • Smoothly tapered ideal sinc function 461 Wolberg: Image Processing Course NotesInverse Mapping Output image Input image v y u x • Visit output in scanline order • Supersampling approximates area sampling • Popular in computer graphics 462 Wolberg: Image Processing Course NotesForward Mapping Input image Output (accumulator) image v y u x • Visit input in scanline order • Use output accumulator array • 2D antialiasing is difficult • Separable transforms facilitate efficient solution 463 Wolberg: Image Processing Course NotesSeparable Transforms X (u,v),Y (u,v) F(u,v) G(x,v) • F(u, v) is a rowpreserving transformation that maps all input points to their final column positions, i.e., x, v. • G(x, v) is a columnpreserving transformation that maps the x, v points to their final row positions, i.e., x, y. 464 Wolberg: Image Processing Course NotesCatmullSmith Algorithm • First pass Maps image S(u,v) into intermediate image I(x,v) I(x,v) = S(X(u,v), v) • Second pass Maps I(x,v) into target image T(x,y) T(x,y) = I(x, Y(H (v), v)) x where H is the solution to x=X(u,v) for u x 465 Wolberg: Image Processing Course Notes2Pass Rotation F G X, Y 466 Wolberg: Image Processing Course Notes2Pass Perspective F G X, Y 467 Wolberg: Image Processing Course NotesFant’s Algorithm • Forward mapping intensity resampling • Scanline order in input and output • Amenable to hardware implementation 468 Wolberg: Image Processing Course NotesFant’s algorithm: Example (1) XLUT .6 2.3 3.2 3.3 3.9 YLUT 100 106 115 120 I 100 106 92 90 YLUT 100 101 105 113 x I 40 101 106 82 x 469 Wolberg: Image Processing Course NotesFant’s algorithm: Example (2) I (0) (100)((.4)) 40 x .4 .4   I (1) (100) 1 (106) ((1)) 101 x   1.7 1.7   1.4 1.4   I (2) (100) 1 (106) ((.3)) (106)((.7)) 106 x   1.7 1.7   .7 .7   I (3) (106) 1 (92) ((.2)) (92)((.1)) (90)((.6)) 82 x   .9 .9   470 Wolberg: Image Processing Course NotesBibliography • Catmull, E. and A.R. Smith, “3D Transformations of Images in Scanline Order,” Proc. Siggraph ‘80, pp. 279285, 1980. • Fant, Karl M., “A Nonaliasing, RealTime Spatial Transform Technique,” IEEE Computer Graphics and Applications, vol. 6, no. 3, pp. 7180, 1986. • Wolberg, George, Digital Image Warping, IEEE Computer Society Press, Los Alamitos, CA 1990. 471 Wolberg: Image Processing Course NotesImage Morphing Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review digital image morphing, a powerful visual effect for creating a fluid transformation from one image to another. Background Morphing Algorithms • Mesh (parametric grid) • Line pairs • Scattered points Examples 473 Wolberg: Image Processing Course NotesInbetween Image Generation Crossdissolve vs. Image Morphing crossdissolve source dest. image image warped warped source dest. image image image morphing 474 Wolberg: Image Processing Course NotesMesh Warp (1) 475 Wolberg: Image Processing Course NotesMesh Warp (2) Source / Target Meshes 476 Wolberg: Image Processing Course NotesMesh Warp (3) Intermediate Mesh 477 Wolberg: Image Processing Course NotesMesh Warp (4) Horizontal Resampling 478 Wolberg: Image Processing Course NotesMesh Warp (5) Vertical Resampling 479 Wolberg: Image Processing Course NotesLine Warp (1) Q’ Q v v X X’ u u PP’ Destination Image Source Image (X P) Perp(Q P) (X P)(Q P) v u 2 Q P Q P v Perp(Q P) X P u(Q P) Q P 480 Wolberg: Image Processing Course NotesLine Warp (2) ’ Q 1 Q 1 Q’ 2 Q 2 X v v 1 2 v 1 v 2 X u 2 X’ 1 u 2 u uX’ 1 1 X’ 2 ’ P 2 P 2 ’ P P 1 1 Destination Image Source Image 481 Wolberg: Image Processing Course NotesScattered Point Constraints Treat correspondence data as scattered data: x,y = X(u,v), Y(u,v) X(u,v) v u 482 Wolberg: Image Processing Course NotesExample 483 Wolberg: Image Processing Course NotesUniform Transition 484 Wolberg: Image Processing Course NotesNonuniform Transition 485 Wolberg: Image Processing Course NotesMorph Sequence 486 Wolberg: Image Processing Course NotesMorph Sequence 487 Wolberg: Image Processing Course NotesBibliography • Beier, T. and S. Neely, “FeatureBased Image Metamorphosis,” Proc. Siggraph ‘92, pp. 3542, 1992. • Smythe, D., “A TwoPass Mesh Warping Algorithm for Object Transformation and Image Interpolation,” ILM Technical Memo 1030, Computer Graphics Dept., Lucasfilm Ltd., 1990. • Lee, S., K.Y. Chwa, S.Y. Shin, and G. Wolberg, “Image Metamorphosis Using Snakes and FreeForm Deformations,” Proc. Siggraph ‘95, pp. 439448, 1992. • Lee, S., G. Wolberg, K.Y. Chwa, and S.Y. Shin, “Image Metamorphosis with Scattered Feature Constraints,” IEEE Trans. Visualization and Computer Graphics, vol. 2, no. 4, pp. 337 354, 1996. • Wolberg, G., “Image Morphing: A Survey”, Visual Computer, vol. 14, no. 8/9, pp. 360372, 1998. • Wolberg, G., Digital Image Warping, IEEE Computer Society Press, Los Alamitos, CA 1990. 488 Wolberg: Image Processing Course NotesRobust LogPolar Registration Prof. George Wolberg Dept. of Computer Science City College of New YorkObjectives • In this lecture we review image registration, a process to align multiple images together. Affine registration Perspective registration Polar coordinates Logpolar algorithm Examples 490 Wolberg: Image Processing Course NotesIntroduction • Image registration refers to the problem of aligning a set of images. • The images may be taken at different times, by different sensors, or from different viewpoints. 491 Wolberg: Image Processing Course NotesApplications • Multisensor data fusion integrating information taken from different sensors • Image analysis / surveillance / change detection for images taken at different times/conditions • Image mosaics generating large panoramas from overlapping images • Superresolution images / noise removal integrating multiple registered images 492 Wolberg: Image Processing Course NotesGeometric Transformations Geometric transformation models for registration: • Rigid transformation translation, rotation (3 parameters) • Affine transformation translation, rotation, scale, shear (6 parameters) • Perspective transformation affine, perspective foreshortening (8 parameters) • Local transformation terrain relief, nonlinear, nonrigid 493 Wolberg: Image Processing Course NotesParameter Estimation (1) Moderate perspective Severe perspective o o o β, γ 40 40 β, γ 90 494 Wolberg: Image Processing Course NotesParameter Estimation (2) The objective function (similarity measure): N 2 ' ' 2  (a) (I (x, y) I (Q x , y ))  1 2 a i1 Q is the affine or perspective transform applied to image I We use the LevenbergMarquardt method to 2. minimize the above objective function for each pyramid level. 495 Wolberg: Image Processing Course NotesAffine Registration Image I Image I I registered to I 1 2 2 1 ° Actual parameters: RST = (30 , 1.5, 30, 30) ° Estimated parameters: RST = (29.99 , 1.51, 29.51, 30.66) 496 Wolberg: Image Processing Course NotesLocal Affine Approximation 2 X 2 tiles 4 X 4 tiles 497 Wolberg: Image Processing Course NotesPiecewise Affine Approximation • Consider I subdivided into a regular grid of tiles 2 i i • For each tile T search for similar tile T in I by 2 1 1 performing affine registration using LMA • Get collection of affine parameters and tile centers • Plug into system of equations • Solve for 8 persp. parameters using the pseudoinverse 498 Wolberg: Image Processing Course NotesExample (1) 499 Wolberg: Image Processing Course NotesExample (1) 500 Wolberg: Image Processing Course NotesPerspective Registration 501 Wolberg: Image Processing Course NotesGlare Removal • The camera's flash introduces a highlight. • Each highlight will appear at a different location on each image when captured from different viewpoints. • If we register these images onto a common frame, we can render the final image by fusing the registered images. 502 Wolberg: Image Processing Course NotesExample (2) •Video mosaic 503 Wolberg: Image Processing Course NotesExample (3) •Video stabilization 504 Wolberg: Image Processing Course NotesDiscussion • Method: LevenbergMarquadt Alg. (nonlinear least squares) • Benefit: estimation of realvalued parameters • Drawback: the registered images must be fairly close in ° scale (1.5), rotation (45 ), and translation • Solution: logpolar registration 505 Wolberg: Image Processing Course NotesPolar Coordinates Radial lines in (x,y) Cartesian space map to horizontal lines in (r,) polar coordinate space. 2 2 r (x x ) (y y ) c c y y 1 c  tan x x c 506 Wolberg: Image Processing Course NotesPolar Coordinates: Example Input image Polar transformation Output in Cartesian space Circular shift in polar space 507 Wolberg: Image Processing Course NotesLogPolar Coordinates Radial lines in (x,y) map to horizontal lines in (log r, ) 508 Wolberg: Image Processing Course NotesExample Image I Logpolar(I ) 1 1 ° Image I : 3x; 45 Logpolar(I ) 2 2 509 Wolberg: Image Processing Course NotesLogPolar Registration • Benefit: Correlation recovers rotation/scale (fast) • Drawback: Image origins must be known • Solution: ' 1. Crop central region I from I 1 1 ' ' 2. Compute I , the log polar transformation of I 1p 1 3. For all positions (x, y) in I : 2 ' Crop region I 2 ' Compute I 2 p ' ' Cross correlateI and I (dx, dy) 1p 2 p if maximum correlation, save (x, y) and (dx, dy) 4. Scale dx 5. Rotation dy 6. Translation (x, y) 510 Wolberg: Image Processing Course NotesScanning Window 511 Wolberg: Image Processing Course NotesExample 512 Wolberg: Image Processing Course NotesExample  Rotation = 21 Scale = 0.469 Translation = (11,51) 513 Wolberg: Image Processing Course NotesBiological Motivations (1) • Despite the difficulties of nonlinear processing, the log polar transform has received considerable attention. Marshal41 et. al. and Hubel74 discovered a log polar mapping in the primate visual system and this conformal map is an accepted model of the representation of the retina in the primary visual cortex in primates Schwartz79 Weinshall87. 514 Wolberg: Image Processing Course NotesBiological Motivations (2) The left figure shows contours of cone cell density. In fact, the density of the ganglion cells which transmit information out of the retina is distributed logarithmically. From Osterberg, G. (1935) Benefits: 1. Reduce the amount of information traversing the optical nerve while maintaining high resolution in the fovea and capturing a wide field of view in the periphery. 2. Invariant to scale and rotation. 515 Wolberg: Image Processing Course NotesBiological Motivations (3) Fovea based vision 516 Wolberg: Image Processing Course NotesRelated Work • The FourierMellin transform uses logpolar mapping to align images related by scale, rotation, and translation Casasent76, Decastro87, Chen94, Reddy96, Lucchese96, Chang97, Lucchese97a, Lucchese97b, Stone99, Stone03, Keller03. • Detect lines in logHough space Weiman79, Giesler98, Weiman90, Young00. • Recognizing 2D objects that were arbitrarily rotated or scaled Sandini92, Ferrari95. • Tracking a moving object Capurro97. • Estimation of timetoimpact from optical flow Sandini91 • Finding disparity map in stereo images Sandini01. • Using logpolar transform to calculate timetocrash for mobile vehicles Pardo02. • A foveated binocular stereo system using log polar transforms Bernardino02. 517 Wolberg: Image Processing Course NotesPrevious work: FourierMellin I (x, y) I (s(x cos( ) y sin( ) x ), s(xsin( ) y cos( )) y ) 1 2 0 0 0 0 0 0  cos( ) sin( ) sin( ) cos( ) 1 x 0 y 0 x 0 y 0 J ( x y ) 0 x 0 y F ( , ) F ( , )e 1 x y 2 s s s • The magnitude of spectra F is a rotated and 1 scaled replica of F . We can recover this 2 rotation and scale by representing the spectra F and F in logpolar coordinates: 1 2 1 F (logr, ) F (logr log s, ) 1 2 0 s • Phasecorrelation can be used to recover the amount shift 518 Wolberg: Image Processing Course NotesPrevious work: FourierMellin ◦ • Example: s = 1.5, θ = 36 519 Wolberg: Image Processing Course NotesPrevious work: FourierMellin •Benefit: We can search for the scale and rotation independent of the translations. •Drawbacks: Large translation, scale, or mild perspective will alter the coefficients of the finite and discrete Fourier transform. 520 Wolberg: Image Processing Course NotesPrevious work: FourierMellin (a) Reference image (b) Target image (real) (c) Target image (synthetic) Power spectra of (a), (b), and (c) 521 Wolberg: Image Processing Course NotesExample from INRIA 522 Wolberg: Image Processing Course NotesExample (1) 523 Wolberg: Image Processing Course NotesExample (2) 524 Wolberg: Image Processing Course NotesExample (3) Recovered parameters: ◦ θ = 228.17 , S = 3.528217 Reference image Target image Tx=32 Ty=35 525 Wolberg: Image Processing Course NotesExample: Large Perspective (1) 526 Wolberg: Image Processing Course NotesExample: Large Perspective (2) 527 Wolberg: Image Processing Course NotesSummary •Logpolar registration: recovers largescale rotation/scale/translation applies fast correlation in logpolar space exploits fast multiresolution processing initial estimate for affine/perspective registration •Perspective registration: handles small distortions with subpixel accuracy applies LevenbergMarquardt algorithm nonlinear least squares optimization exploits fast multiresolution processing 528 Wolberg: Image Processing Course Notes
sharer
Presentations
Free
Document Information
Category:
Presentations
User Name:
VoiletFord
User Type:
Professional
Country:
United States
Uploaded Date:
12-07-2017