[iadftscaleproperty] [Up] [iahotelling] Lessons

iaconvteo
Illustrate the convolution theorem

Description

In this demonstration, the convolution theorem is illustrated by comparing the result of the convolution with its equivalent filtering in frequency domain.

Demo Script

Reading and ROI selection

The image is read and displayed 

>>> fin = iaread('lenina.pgm')

                  
>>> iashow(fin)

                  
(256, 256) Min= 36 Max= 255 Mean=88.958 Std=46.38

                  
>>> froi = iaroi(fin, (90,70), (200,180))

                  

                  
>>> iashow(froi)

                  
(111, 111) Min= 45 Max= 209 Mean=112.001 Std=29.94

                

                

                  
                    

                      
                      
                      
                    

                    

                      fin
                      froi
                      
                    

                  
                

              

            

          

        

        

          
Convolution with the Laplacian kernel

          

            

              

                        The image is convolved (periodicaly) with the 3x3 Laplacian kernel
                    
              

              

                

                  
>>> import Numeric

                  

                  
>>> fd = froi.astype(Numeric.Float)

                  

                  
>>> h = Numeric.array([[-2,-1,0],[-1,0,1],[0,1,2]])

                  

                  
>>> g = iapconv(fd,h)

                  

                  
>>> iashow(g)

                  
(111, 111) Min= -308.0 Max= 391.0 Mean=0.000 Std=41.56

                

                

                  
                    

                      
                      
                    

                    

                      g
                      
                    

                  
                

              

            

          

        

        

          
Equivalent filter in frequency domain

          

            

              

                        The 3x3 kernel is zero padded to the size of the input image and periodicaly translated so that the center of the kernel stays at the top-left image corner. Its spectrum is visualized.
                    
              

              

                

                  
>>> hx = Numeric.zeros(Numeric.array(froi.shape))

                  

                  
>>> hx[:h.shape[0],:h.shape[1]] = h

                  

                  
>>> hx = iaptrans(hx,-Numeric.floor((Numeric.array(h.shape)-1)/2).astype(Numeric.Int))

                  

                  
>>> H = iadft(hx);

                  

                  
>>> iashow(iadftview(H))

                  
(111, 111) Min= 0 Max= 254 Mean=147.832 Std=66.87

                

                

                  
                    

                      
                      
                    

                    

                      iadftview(H)
                      
                    

                  
                

              

            

          

        

        

          
Filtering in the frequency domain

          

            

              

                        The image is filtered by multiplying its DFT by the frequency mask computed in the previous step.
                    
              

              

                

                  
>>> F = iadft(fd)

                  

                  
>>> G = F * H

                  

                  
>>> print "Is symmetrical:", iaisdftsym(G)

                  
Is symmetrical: 1

                  
>>> iashow(iadftview(G))

                  
(111, 111) Min= 0 Max= 255 Mean=144.841 Std=41.62

                  
>>> g_aux = iaidft(G).real

                  

                  
>>> iashow(g_aux)

                  
(111, 111) Min= -308.0 Max= 391.0 Mean=-0.000 Std=41.56

                

                

                  
                    

                      
                      
                      
                    

                    

                      iadftview(G)
                      g_aux
                      
                    

                  
                

              

            

          

        

        

          
Comparing the results

          

            

              

                        Both images, filtered by the convolution and filtered in the frequency domain are compared to see that they are the same. The small differencies are due to numerical precision errors.
                    
              

              

                

                  
>>> e = abs(g - g_aux)

                  

                  
>>> print "Max error:", max(Numeric.ravel(e))

                  
Max error: 4.15525391873e-011

                

              

            

          

        

      

    

    

      

        

          
                    [iadftscaleproperty]
                
                    [Up]
                    
                    [iahotelling]
                
          http://www.python.org