数字图像处理-冈萨雷斯-课件英文版-Chapter05-图像复原.ppt

1、Digital Image ProcessingChapter 5:Image Restoration23 June 2006(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image restoration is to restore a degraded image back tothe original image while image enhancement is to manipulate the image so that it is suitable f

2、or a specificapplication.Degradation model:),(),(),(),(yxyxhyxfyxgwhere h(x,y)is a system that causes image distortion and(x,y)is noise.Noise cannot be predicted but can be approximately described instatistical way using the probability density function(PDF)Gaussian noise:222/)(21)(zezpRayleigh nois

3、eazfor 0for )(2)(/)(2azeazbzpbazErlang(Gamma)noise0zfor 00for )()!1()(1zeazbzazpazbbExponential noiseUniform noiseImpulse(salt&pepper)noiseazaezp)(otherwise 0afor a-b1)(bzzpotherwise 0for for )(bzPazPzpba(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.PDF tells

4、 how mucheach z value occurs.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Original imageHistogramDegraded images),(),(),(yxyxfyxg(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Original imageHistogramDegraded images),()

5、,(),(yxyxfyxg(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Periodic noise looks like dotsIn the frequencydomain(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.We cannot use the image histogram to estimate noise PDF.It is

6、 better to use the histogram of one areaof an image that has constant intensity to estimate noise PDF.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Band reject filterRestored imageDegraded imageDFTPeriodic noisecan be reduced bysetting frequencycomponentscorr

7、esponding to noise to zero.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Use to eliminate frequency components in some bandsPeriodic noise from theprevious slide that is Filtered out.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2

8、nd Edition.A notch reject filter is used to eliminate some frequency components.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Degraded imageDFTNotch filter(freq.Domain)Restored imageNoise(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processi

9、ng,2nd Edition.Degraded imageDFT(no shift)Restored imageNoiseDFT of noiseArithmetic mean filter or moving average filter(from Chapter 3)xyStstsgmnyxf),(),(1),(Geometric mean filtermnStsxytsgyxf1),(),(),(mn=size of moving windowDegradation model:),(),(),(),(yxyxhyxfyxgTo remove this part(Images from

10、Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Original imageImage corrupted by AWGNImage obtained using a 3x3geometric mean filterImage obtained using a 3x3arithmetic mean filterAWGN:Additive White Gaussian NoiseHarmonic mean filterxyStstsgmnyxf),(),(1),(Contraharmonic me

11、an filterxyxyStsQStsQtsgtsgyxf),(),(1),(),(),(mn=size of moving windowWorks well for salt noisebut fails for pepper noiseQ=the filter orderPositive Q is suitable for eliminating pepper noise.Negative Q is suitable for eliminating salt noise.For Q=0,the filter reduces to an arithmetic mean filter.For

12、 Q=-1,the filter reduces to a harmonic mean filter.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by pepper noise with prob.=0.1Image corrupted by salt noise with prob.=0.1Image obtained using a 3x3contra-harmonic mean filterWith Q=1.5Image obt

13、ained using a 3x3contra-harmonic mean filterWith Q=-1.5(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by pepper noise with prob.=0.1Image corrupted by salt noise with prob.=0.1Image obtained using a 3x3contra-harmonic mean filterWith Q=-1.5Imag

14、e obtained using a 3x3contra-harmonic mean filterWith Q=1.5subimageOriginal imageMoving windowStatistic parametersMean,Median,Mode,Min,Max,Etc.Output imageMedian filter),(median),(),(tsgyxfxyStsMax filter),(max),(),(tsgyxfxyStsMin filter),(min),(),(tsgyxfxyStsMidpoint filter),(min),(max21),(),(),(ts

15、gtsgyxfxyxyStsStsReduce“dark”noise (pepper noise)Reduce“bright”noise (salt noise)A median filter is good for removing impulse,isolated noiseDegraded imageSalt noisePepper noiseMovingwindowSorted arraySalt noisePepper noiseMedianFilter outputNormally,impulse noise has high magnitude and is isolated.W

16、hen we sort pixels in the moving window,noise pixels are usually at the ends of the array.Therefore,its rare that the noise pixel will be a median value.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by salt-and-pepper noise with pa=pb=0.1Image

17、s obtained using a 3x3 median filter1423(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by pepper noise with prob.=0.1Image corrupted by salt noise with prob.=0.1Image obtained using a 3x3max filterImage obtained using a 3x3min filterxyStsrtsgdm

18、nyxf),(),(1),(where gr(s,t)represent the remaining mn-d pixels after removing the d/2 highest and d/2 lowest values of g(s,t).This filter is useful in situations involving multiple typesof noise such as a combination of salt-and-pepper and Gaussian noise.Formula:(Images from Rafael C.Gonzalez and Ri

19、chard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by additiveuniform noiseImage obtained using a 5x5arithmetic mean filterImage additionallycorrupted by additivesalt-and-pepper noise122Image obtained using a 5x5geometric mean filter2Image corrupted by additiveuniform noiseImage obtai

20、ned using a 5x5 median filterImage additionallycorrupted by additivesalt-and-pepper noise122Image obtained using a 5x5alpha-trimmed mean filterwith d=52Image obtained using a 5x5arithmetic mean filterImage obtained using a 5x5geometric mean filterImage obtained using a 5x5 median filterImage obtaine

21、d using a 5x5alpha-trimmed mean filterwith d=5-Filter behavior depends on statistical characteristics of local areas inside mxn moving window-More complex but superior performance compared with“fixed”filtersStatistical characteristics:General concept:xyStsLtsgmnm),(),(1Local mean:Local variance:xySt

22、sLLmtsgmn),(22),(12Noise variance:Purpose:want to preserve edges1.If 2 is zero,No noisethe filter should return g(x,y)because g(x,y)=f(x,y)2.If L2 is high relative to 2,Edges(should be preserved),the filter should return the value close to g(x,y)3.If L2=2,Areas inside objectsthe filter should return

23、 the arithmetic mean value mLLLmyxgyxgyxf),(),(),(22Formula:Concept:(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by additiveGaussian noise with zero meanand 2=1000Imageobtained using a 7x7arithmeticmean filterImageobtained using a 7x7geometri

24、cmean filterImageobtained using a 7x7adaptivenoise reduction filterAlgorithm:Level A:A1=zmedian zminA2=zmedian zmaxIf A1 0 and A2 0,goto level BElse increase window sizeIf window size 0 and B2 0,return zxyElse return zmedianzmin=minimum gray level value in Sxyzmax=maximum gray level value in Sxyzmed

25、ian=median of gray levels in Sxyzxy=gray level value at pixel(x,y)Smax=maximum allowed size of SxywherePurpose:want to remove impulse noise while preserving edgesLevel A:A1=zmedian zmin A2=zmedian zmax Else Window is not big enough increase window sizeIf window size 0 and B2 0 and A2 0,goto level BL

26、evel B:Determine whether zmedianis an impulse or notDetermine whether zxyis an impulse or not(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Image corrupted by salt-and-pepper noise with pa=pb=0.25Image obtained using a 7x7median filterImage obtained using an a

27、daptivemedian filter withSmax=7More small details are preserved(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Degradation model:),(),(),(),(yxyxhyxfyxgPurpose:to estimate h(x,y)or H(u,v),(),(),(),(vuNvuHvuFvuGMethods:1.Estimation by Image Observation2.Estimati

28、on by Experiment3.Estimation by ModelingorWhy?If we know exactly h(x,y),regardless of noise,we can do deconvolution to get f(x,y)back from g(x,y).f(x,y)f(x,y)*h(x,y)g(x,y)SubimageReconstructedSubimage),(vuGs),(yxgs),(yxfsDFTDFT),(vuFsRestoration process byestimationOriginal image(unknown)Degraded im

29、age),(),(),(),(vuFvuGvuHvuHsssEstimated Transfer functionObservationThis case is used when weknow only g(x,y)and cannot repeat the experiment!Used when we have the same equipment set up and can repeat the experiment.Input impulse imageSystemH()Response image fromthe system),(vuG),(yxg),(yxAAyxADFT),

30、(AvuGvuH),(),(DFTDFT(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Used when we know physical mechanism underlying the image formation process that can be expressed mathemat

31、ically.AtmosphericTurbulence model6/522)(),(vukevuHExample:Original imageSevere turbulencek=0.00025k=0.001k=0.0025Low turbulenceMild turbulenceAssume that camera velocity is)(),(00tytxThe blurred image is obtained bydttyytxxfyxgT)(),(),(000where T=exposure time.dtdxdyetyytxxfdxdyedttyytxxfdxdyeyxgvu

32、GTvyuxjvyuxjTvyuxj 0)(200)(2000)(2)(),()(),(),(),(dtevuFdtevuFdtdxdyetyytxxfvuGTtvytuxjTtvytuxjTvyuxj 0)()(20)()(20)(2000000),(),()(),(),(Then we get,the motion blurring transfer function:dtevuHTtvytuxj0)()(200),(For constant motion),()(),(00btattytx)(0)(2)(sin()(),(vbuajTvbuajevbuavbuaTdtevuH(Image

33、s from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.For constant motion)()(sin()(),(vbuajevbuavbuaTvuHOriginal imageMotion blurred imagea=b=0.1,T=1(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.after we obtain H(u,v),we can estimat

34、e F(u,v)by the inverse filter:),(),(),(),(),(),(vuHvuNvuFvuHvuGvuFFrom degradation model:),(),(),(),(vuNvuHvuFvuGNoise is enhancedwhen H(u,v)is small.To avoid the side effect of enhancing noise,we can apply this formulation to ponent(u,v)with in a radius D0 from the center of H(u,v).In practical,the

35、 inverse filter is notPopularly used.6/522)(0025.0),(vuevuHOriginal imageBlurred imageDue to TurbulenceResult of applyingthe full filterResult of applyingthe filter with D0=70 Result of applyingthe filter with D0=40 Result of applyingthe filter with D0=85(Images from Rafael C.Gonzalez and Richard E.

36、Wood,Digital Image Processing,2nd Edition.Objective:optimize mean square error:22)(ffEe),(),(/),(),(),(),(1 ),(),(/),(),(),(),(),(),(),(),(),(),(222*2*vuGvuSvuSvuHvuHvuHvuGvuSvuSvuHvuHvuGvuSvuHvuSvuSvuHvuFffffWiener Filter Formula:whereH(u,v)=Degradation functionS(u,v)=Power spectrum of noiseSf(u,v)

37、=Power spectrum of the undegraded image),(),(/),(),(),(),(1),(22vuGvuSvuSvuHvuHvuHvuFfWiener Filter Formula:Approximated Formula:),(),(),(),(1),(22vuGKvuHvuHvuHvuFDifficult to estimatePractically,K is chosen manually to obtained the best visual result!Original imageBlurred imageDue to TurbulenceResu

38、lt of the full inverse filterResult of the inverse filter with D0=70 Result of the full Wiener filterOriginal imageResult of the inverse filter with D0=70 Result of the Wiener filterBlurred imageDue to Turbulence(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.I

39、mage degradedby motion blur+AWGNResult of theinverse filterResult of theWiener filter2=6502=3252=130Note:K is chosenmanuallyDegradation model:),(),(),(),(yxyxhyxfyxgWritten in a matrix formHfgObjective:to find the minimum of a criterion function101022),(MxNyyxfCSubject to the constraint22fHg),(),(),

40、(),(),(22*vuGvuPvuHvuHvuFWe get a constrained least square filterwhereP(u,v)=Fourier transform of p(x,y)=010141010wherewwwT2),(),(),(),(),(22*vuGvuPvuHvuHvuFConstrained least square filter is adaptively adjusted to achieve the best result.Results from the previous slide obtained from the constrained

41、 least square filter Image degradedby motion blur+AWGNResult of theConstrainedLeast square filterResult of theWiener filter2=6502=3252=130DefinefHgrIt can be shown that2)(rrrTWe want to adjust gamma so thata22rwhere a=accuracy factor1.Specify an initial value of 2.Compute3.Stop if is satisfiedOtherw

42、ise return step 2 after increasing if or decreasing ifUse the new value of to recompute112ra22ra22r),(),(),(),(),(22*vuGvuPvuHvuHvuF),(),(),(),(),(22*vuGvuPvuHvuHvuF),(),(),(),(vuFvuHvuGvuR101022),(1MxNyyxrMNr101022),(1MxNymyxMN1010),(1MxNyyxMNmmMN2222rFor computingFor computing(Images from Rafael C

43、.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Original imageBlurred imageDue to TurbulenceResults obtained from constrained least square filtersUse wrong noise parametersCorrect parameters:Initial =10-5Correction factor=10-6a=0.252=10-5Wrong noise parameter2=10-2Use correct noise

44、 parameters),(),(),(),(),(),(),(),(12*2*vuGvuSvuSvuHvuHvuHvuHvuFf This filter represents a family of filters combined into a single expression=1 the inverse filter=0 the Parametric Wiener filter=0,=1 the standard Wiener filter=1,0.5 More like the Wiener filterAnother name:the spectrum equalization f

45、ilterThese transformations are often called rubber-sheet transformations:Printing an image on a rubber sheet and then stretch this sheet accordingto some predefine set of rules.A geometric transformation consists of 2 basic operations:1.A spatial transformation:Define how pixels are to be rearranged

46、 in the spatiallytransformed image.2.Gray level interpolation:Assign gray level values to pixels in the spatiallytransformed image.Distorted image g1.Select coordinate(x,y)in f to be restored2.Compute),(yxrx),(yxsy 3.Go to pixel in a distorted image g),(yxImage f to be restored4.get pixel value atBy

47、 gray level interpolation),(yxg5.store that value in pixel f(x,y)135),(yx),(yx(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.To map between pixel coordinate(x,y)of f and pixel coordinate(x,y)of g),(yxrx),(yxsy For a bilinear transformation mapping between a pa

48、ir of Quadrilateral regions4321),(cxycycxcyxrx8765),(cxycycxcyxsy),(yx),(yxTo obtain r(x,y)and s(x,y),we need to know 4 pairs of coordinates and its correspondingwhich are called tiepoints.),(yx),(yx(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Since may not

49、be at an integer coordinate,we need to Interpolate the value of ),(yx),(yxgExample interpolation methods that can be used:1.Nearest neighbor selection2.Bilinear interpolation3.Bicubic interpolationOriginal image and tiepointsTiepoints of distortedimageDistorted imageRestored imageUse nearestneighbor

50、 intepolation(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.Original image and tiepointsTiepoints of distortedimageDistorted imageRestored imageUse bilinear intepolation(Images from Rafael C.Gonzalez and Richard E.Wood,Digital Image Processing,2nd Edition.(Ima