In [1]:
import math
import os
import rasterio
from sklearn.preprocessing import *
from skimage.measure import *
from __future__ import division
import cv2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib qt

font = {'family': 'serif',
        'color':  'black',
        'weight': 'normal',
        'size': 12,
        }

In [2]:
snow_compare = ['/home/cparr/water_tracks/sub_wt_2012.tif',
                '/home/cparr/water_tracks/sub_wt_2013.tif',
                '/home/cparr/water_tracks/sub_wt_2015.tif']

src_dict_2012 = {}
src_dict_2013 = {}
src_dict_2015 = {}

dcts = [src_dict_2012, src_dict_2013, src_dict_2015]

i = 0
for f in snow_compare:

    dcts[i]['name'] = os.path.basename( f ).replace( '.tif', '_src' )
    src = rasterio.open( f )
    meta = src.meta
    
    dcts[i]['actual snow depth'] = np.ma.masked_values( src.read(1), src.nodata )
    dcts[i]['mean snow depth'] = np.mean( dcts[i]['actual snow depth'] )
    dcts[i]['std snow depth'] = np.std( dcts[i]['actual snow depth'] )
    
    dcts[i]['blur snow depth'] = cv2.GaussianBlur(dcts[i]['actual snow depth'], (3,3), 0)
    dcts[i]['mean blur depth'] = np.mean( dcts[i]['blur snow depth'] )
    dcts[i]['std blur depth'] = np.std( dcts[i]['blur snow depth'] )
    
    dcts[i]['norm blur snow depth'] = maxabs_scale(dcts[i]['blur snow depth'], axis=1, copy=True)
    dcts[i]['mean norm blur depth'] = np.mean( dcts[i]['norm blur snow depth'] )
    dcts[i]['std norm blur depth'] = np.std( dcts[i]['norm blur snow depth'] )
        
    dcts[i]['norm drift mask'] = np.ma.MaskedArray(dcts[i]['norm blur snow depth'] > 0.5)
    dcts[i]['norm drift percent'] = ( np.sum( dcts[i]['norm drift mask'] ) / np.size( dcts[i]['norm drift mask'] )) * 100
    
    dcts[i]['norm scour mask'] = np.ma.MaskedArray(dcts[i]['norm blur snow depth'] <= 0.5)
        
    dcts[i]['norm drift snow'] = dcts[i]['norm drift mask'] * dcts[i]['norm blur snow depth']
    dcts[i]['norm drift snow'][dcts[i]['norm drift snow'] == 0] = np.nan
    dcts[i]['mean norm drift depth'] = np.nanmean( dcts[i]['norm drift snow'] )
    dcts[i]['std norm drift depth'] = np.nanstd( dcts[i]['norm drift snow'] )
    
    dcts[i]['actual drift snow'] = dcts[i]['norm drift mask'] * dcts[i]['actual snow depth']
    dcts[i]['actual drift snow'][dcts[i]['actual drift snow'] == 0] = np.nan
    dcts[i]['mean actual drift depth'] = np.nanmean( dcts[i]['actual drift snow'] )
    dcts[i]['std actual drift depth'] = np.nanstd( dcts[i]['actual drift snow'] )
    
    dcts[i]['norm scoured snow'] = dcts[i]['norm scour mask'] * dcts[i]['norm blur snow depth']
    dcts[i]['actual scour snow'] = dcts[i]['norm scour mask'] * dcts[i]['actual snow depth']
        
    i += 1

In [ ]:
# Actual Depths

plt.suptitle( 'lidar snow depths [m]', x = 0.75, y = 0.33, fontsize = 20, backgroundcolor = 'powderblue' )

plt.subplot( 2,2,1 )
plt.imshow( src_dict_2012['actual snow depth'] )
plt.title( "2012: mean " + str( src_dict_2012['mean snow depth'] )[0:4] + 
          ", std " + str( src_dict_2012['std snow depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,2 )
plt.imshow( src_dict_2013['actual snow depth'] )
plt.title( "2013: mean " + str( src_dict_2013['mean snow depth'] )[0:4] + 
          ", std " + str( src_dict_2013['std snow depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,3 )
plt.imshow( src_dict_2015['actual snow depth'] )
plt.title( "2015: mean " + str( src_dict_2015['mean snow depth'] )[0:4] + 
          ", std " + str( src_dict_2015['std snow depth'] )[0:4], font )
plt.colorbar()

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/mean_snow.png', dpi = 300)

In [ ]:
# Gaussian Blur

plt.suptitle( '3x3 Gaussian Blur [m]', x = 0.75, y = 0.33, fontsize = 18, backgroundcolor = 'powderblue' )

plt.subplot( 2,2,1 )
plt.imshow( src_dict_2012['blur snow depth'] )
plt.title( "2012: mean " + str( src_dict_2012['mean blur depth'] )[0:4] + 
          ", std " + str( src_dict_2012['std blur depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,2 )
plt.imshow( src_dict_2013['blur snow depth'] )
plt.title( "2013: mean " + str( src_dict_2013['mean blur depth'] )[0:4] + 
          ", std " + str( src_dict_2013['std blur depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,3 )
plt.imshow( src_dict_2015['blur snow depth'] )
plt.title( "2015: mean " + str( src_dict_2015['mean blur depth'] )[0:4] + 
          ", std " + str( src_dict_2015['std blur depth'] )[0:4], font )
plt.colorbar()

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/blur_snow.png', dpi = 300)

In [38]:
# Normal Blur

plt.suptitle( 'Normalized Blur', x = 0.75, y = 0.33, fontsize = 18, backgroundcolor = 'powderblue' )

plt.subplot( 2,2,1 )
plt.imshow( src_dict_2012['norm blur snow depth'] )
plt.title( "2012: mean " + str( src_dict_2012['mean norm blur depth'] )[0:4] + 
          ", std " + str( src_dict_2012['std norm blur depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,2 )
plt.imshow( src_dict_2013['norm blur snow depth'] )
plt.title( "2013: mean " + str( src_dict_2013['mean norm blur depth'] )[0:4] + 
          ", std " + str( src_dict_2013['std norm blur depth'] )[0:4], font )
plt.colorbar()

plt.subplot( 2,2,3 )
plt.imshow( src_dict_2015['norm blur snow depth'] )
plt.title( "2015: mean " + str( src_dict_2015['mean norm blur depth'] )[0:4] + 
          ", std " + str( src_dict_2015['std norm blur depth'] )[0:4], font )
plt.colorbar()

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/norm_blur_snow.png', dpi = 300)

In [ ]:
# Drift Mask %

plt.suptitle( 'Drift Mask (Norm > 0.5)', x = 0.75, y = 0.33, fontsize = 20, backgroundcolor = 'powderblue' )

plt.subplot( 2,2,1 )
plt.imshow( src_dict_2012['norm drift mask'] , cmap = 'gray' )
plt.title( "2012 Drift % " + str( src_dict_2012['norm drift percent'] )[0:4])
          
plt.subplot( 2,2,2 )
plt.imshow( src_dict_2013['norm drift mask'] , cmap = 'gray' )
plt.title( "2013 Drift % " + str( src_dict_2013['norm drift percent'] )[0:4])

plt.subplot( 2,2,3 )
plt.imshow( src_dict_2015['norm drift mask'] , cmap = 'gray' )
plt.title( "2015 Drift % " + str( src_dict_2015['norm drift percent'] )[0:4])

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/drift_mask.png', dpi = 300)

In [ ]:
# Drift Snow

font = {'family': 'serif',
        'color':  'black',
        'weight': 'normal',
        'size': 10,
        }

plt.subplot( 2,3,1 )
plt.imshow( src_dict_2012['norm drift snow'] )
plt.title( "2012: mean " + str( src_dict_2012['mean norm drift depth'] )[0:4] + 
          ", std " + str( src_dict_2012['std norm drift depth'] )[0:4], font )
          
plt.subplot( 2,3,2 )
plt.imshow( src_dict_2013['norm drift snow'] )
plt.title( "2013: mean " + str( src_dict_2013['mean norm drift depth'] )[0:4] + 
          ", std " + str( src_dict_2013['std norm drift depth'] )[0:4], font )

plt.annotate('Normalized Depths (top) and Actual Depths [m] (bottom)', (0,0), (-75, -33), 
             xycoords='axes fraction', textcoords='offset points', va='top', backgroundcolor = 'powderblue' )

plt.subplot( 2,3,3 )
plt.imshow( src_dict_2015['norm drift snow'] )
plt.title( "2015: mean " + str( src_dict_2015['mean norm drift depth'] )[0:4] + 
          ", std " + str( src_dict_2015['std norm drift depth'] )[0:4], font )

plt.subplot( 2,3,4 )
plt.imshow( src_dict_2012['actual drift snow'] )
plt.title( "2012: mean " + str( src_dict_2012['mean actual drift depth'] )[0:4] + 
          ", std " + str( src_dict_2012['std actual drift depth'] )[0:4], font )
          
plt.subplot( 2,3,5 )
plt.imshow( src_dict_2013['actual drift snow'] )
plt.title( "2013: mean " + str( src_dict_2013['mean actual drift depth'] )[0:4] + 
          ", std " + str( src_dict_2013['std actual drift depth'] )[0:4], font )

plt.subplot( 2,3,6 )
plt.imshow( src_dict_2015['actual drift snow'] )
plt.title( "2015: mean " + str( src_dict_2015['mean actual drift depth'] )[0:4] + 
          ", std " + str( src_dict_2015['std actual drift depth'] )[0:4], font )

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/drifted_snow.png', dpi = 300)

In [52]:
# Mean Square Error Comparison

pair_12_13 = {}
pair_12_15 = {}
pair_13_15 = {}

pair_12_13['actual'] = [ src_dict_2012['actual snow depth'], src_dict_2013['actual snow depth'] ]
pair_12_15['actual'] = [ src_dict_2012['actual snow depth'], src_dict_2015['actual snow depth'] ]
pair_13_15['actual'] = [ src_dict_2013['actual snow depth'], src_dict_2015['actual snow depth'] ]

pair_12_13['actual drift'] = [ src_dict_2012['actual drift snow'], src_dict_2013['actual drift snow'] ]
pair_12_15['actual drift'] = [ src_dict_2012['actual drift snow'], src_dict_2015['actual drift snow'] ]
pair_13_15['actual drift'] = [ src_dict_2013['actual drift snow'], src_dict_2015['actual drift snow'] ]

pair_12_13['norm drift'] = [ src_dict_2012['norm drift snow'], src_dict_2013['norm drift snow'] ]
pair_12_15['norm drift'] = [ src_dict_2012['norm drift snow'], src_dict_2015['norm drift snow'] ]
pair_13_15['norm drift'] = [ src_dict_2013['norm drift snow'], src_dict_2015['norm drift snow'] ]

pair_12_13['norm'] = [ src_dict_2012['norm blur snow depth'], src_dict_2013['norm blur snow depth'] ]
pair_12_15['norm'] = [ src_dict_2012['norm blur snow depth'], src_dict_2015['norm blur snow depth'] ]
pair_13_15['norm'] = [ src_dict_2013['norm blur snow depth'], src_dict_2015['norm blur snow depth'] ]

pair_12_13['intersection'] = src_dict_2012['norm drift mask'].__and__( src_dict_2013['norm drift mask'] )
pair_12_15['intersection'] = src_dict_2012['norm drift mask'].__and__( src_dict_2015['norm drift mask'] )
pair_13_15['intersection'] = src_dict_2013['norm drift mask'].__and__( src_dict_2015['norm drift mask'] )

pair_12_13['either or'] = src_dict_2012['norm drift mask'].__xor__( src_dict_2013['norm drift mask'] )
pair_12_15['either or'] = src_dict_2012['norm drift mask'].__xor__( src_dict_2015['norm drift mask'] )
pair_13_15['either or'] = src_dict_2013['norm drift mask'].__xor__( src_dict_2015['norm drift mask'] )

pairs = [ pair_12_13, pair_12_15, pair_13_15 ]

def mse(imageA, imageB):
    a_mse = ((imageA.astype("float") - imageB.astype("float")) ** 2)
    err = np.nansum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])
    return err, a_mse

for p in pairs:
        p['actual_mse error value'] = mse( p['actual'][0], p['actual'][1] )[0]
        p['actual_mse error array'] = mse( p['actual'][0], p['actual'][1] )[1]
        p['actual drift_mse error value'] = mse( p['actual drift'][0], p['actual drift'][1] )[0]
        p['actual drift_mse error array'] = mse( p['actual drift'][0], p['actual drift'][1] )[1]
        p['norm drift_mse error value'] = mse( p['norm drift'][0], p['norm drift'][1] )[0]
        p['norm drift_mse error array'] = mse( p['norm drift'][0], p['norm drift'][1] )[1]
        p['norm_mse error value'] = mse( p['norm'][0], p['norm'][1] )[0]
        p['norm_mse error array'] = mse( p['norm'][0], p['norm'][1] )[1]

In [48]:
# MSE Actual Figs

plt.subplot( 2,3,1 )
plt.imshow( pair_12_13['actual_mse error array'] )
plt.title( "2012 - 2013 MSE " + str( pair_12_13['actual_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,2 )
plt.imshow( pair_12_15['actual_mse error array'] )
plt.title( "2012 - 2015 MSE  " + str( pair_12_15['actual_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.annotate('Actual Depth Mean Square Error', (0,0), (-33, -33), 
             xycoords='axes fraction', textcoords='offset points', va='top', backgroundcolor = 'powderblue' )

plt.subplot( 2,3,3 )
plt.imshow( pair_13_15['actual_mse error array'] )
plt.title( "2013 - 2015 MSE  " + str( pair_13_15['actual_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

#
          
plt.subplot( 2,3,4 )
plt.imshow( pair_12_13['actual drift_mse error array'] )
plt.title( "2012 - 2013 MSE " + str( pair_12_13['actual drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,5 )
plt.imshow( pair_12_15['actual drift_mse error array'] )
plt.title( "2012 - 2015 MSE  " + str( pair_12_15['actual drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,6 )
plt.imshow( pair_13_15['actual drift_mse error array'] )
plt.title( "2013 - 2015 MSE  " + str( pair_13_15['actual drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/mse_actual.png', dpi = 300)

In [49]:
# MSE norm Figs
     
font = {'family': 'serif',
        'color':  'black',
        'weight': 'normal',
        'size': 8,
        }

plt.subplot( 2,3,1 )
plt.imshow( pair_12_13['norm_mse error array'] )
plt.title( "2012 - 2013 MSE " + str( pair_12_13['norm_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )
          
plt.subplot( 2,3,2 )
plt.imshow( pair_12_15['norm_mse error array'] )
plt.title( "2012 - 2015 MSE  " + str( pair_12_15['norm_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.annotate('Normal Depth Mean Square Error', (0,0), (-33, -33), 
             xycoords='axes fraction', textcoords='offset points', va='top', backgroundcolor = 'powderblue' )

plt.subplot( 2,3,3 )
plt.imshow( pair_13_15['norm_mse error array'] )
plt.title( "2013 - 2015 MSE  " + str( pair_13_15['norm_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )
#     
plt.subplot( 2,3,4 )
plt.imshow( pair_12_13['norm drift_mse error array'] )
plt.title( "2012 - 2013 MSE " + str( pair_12_13['norm drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )
          
plt.subplot( 2,3,5 )
plt.imshow( pair_12_15['norm drift_mse error array'] )
plt.title( "2012 - 2015 MSE  " + str( pair_12_15['norm drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,6 )
plt.imshow( pair_13_15['norm drift_mse error array'] )
plt.title( "2013 - 2015 MSE  " + str( pair_13_15['norm drift_mse error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/mse_norm.png', dpi = 300)

In [15]:
# SSIM

def structural_sim(imageA, imageB):
    struct_sim_result = compare_ssim( imageA, imageB, full = True )
    ssim_value = struct_sim_result[0]
    ssim_arr = struct_sim_result[1]
    return ssim_value, ssim_arr

for p in pairs:
        p['actual_ssim error value'] = structural_sim( p['actual'][0], p['actual'][1] )[0]
        p['actual_ssim error array'] = structural_sim( p['actual'][0], p['actual'][1] )[1]
        p['actual drift_ssim error value'] = structural_sim( p['actual drift'][0], p['actual drift'][1] )[0]
        p['actual drift_ssim error array'] = structural_sim( p['actual drift'][0], p['actual drift'][1] )[1]
        p['norm drift_ssim error value'] = structural_sim( p['norm drift'][0], p['norm drift'][1] )[0]
        p['norm drift_ssim error array'] = structural_sim( p['norm drift'][0], p['norm drift'][1] )[1]
        p['norm_ssim error value'] = structural_sim( p['norm'][0], p['norm'][1] )[0]
        p['norm_ssim error array'] = structural_sim( p['norm'][0], p['norm'][1] )[1]

In [46]:
# SSIM Figs

plt.subplot( 2,3,1 )
plt.imshow( pair_12_13['actual_ssim error array'] )
plt.title( "2012 - 2013 SSIM " + str( pair_12_13['actual_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,2 )
plt.imshow( pair_12_15['actual_ssim error array'] )
plt.title( "2012 - 2015 SSIM  " + str( pair_12_15['actual_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.annotate('Actual (top) & Normalized (bottom) Structural Similarity Index', (0,0), (-100, -33), 
             xycoords='axes fraction', textcoords='offset points', va='top', backgroundcolor = 'powderblue' )

plt.subplot( 2,3,3 )
plt.imshow( pair_13_15['actual_ssim error array'] )
plt.title( "2013 - 2015 SSIM  " + str( pair_13_15['actual_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )          
         
plt.subplot( 2,3,4 )
plt.imshow( pair_12_13['norm_ssim error array'] )
plt.title( "2012 - 2013 SSIM  " + str( pair_12_13['norm_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,5 )
plt.imshow( pair_12_15['norm_ssim error array'] )
plt.title( "2012 - 2015 SSIM  " + str( pair_12_15['norm_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,6 )
plt.imshow( pair_13_15['norm_ssim error array'] )
plt.title( "2013 - 2015 SSIM  " + str( pair_13_15['norm_ssim error value'] )[0:5] )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] ) 

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/ssim.png', dpi = 300)

In [56]:
# Intersection XOR Figs

plt.subplot( 2,3,1 )
plt.imshow( pair_12_13['intersection'] , cmap = 'gray' )
plt.title( "2012 - 2013 Intersection" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,2 )
plt.imshow( pair_12_15['intersection'] , cmap = 'gray' )
plt.title( "2012 - 2015 Intersection" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.annotate('Drift Sets', (0,0), (-0, -33), 
             xycoords='axes fraction', textcoords='offset points', va='top', backgroundcolor = 'powderblue' )

plt.subplot( 2,3,3 )
plt.imshow( pair_13_15['intersection'] , cmap = 'gray' )
plt.title( "2013 - 2015 Intersection" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )          
         
plt.subplot( 2,3,4 )
plt.imshow( pair_12_13['either or'] , cmap = 'gray' )
plt.title( "2012 - 2013 Difference" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,5 )
plt.imshow( pair_12_15['either or'] , cmap = 'gray' )
plt.title( "2012 - 2015 Difference" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] )

plt.subplot( 2,3,6 )
plt.imshow( pair_13_15['either or'] , cmap = 'gray' )
plt.title( "2013 - 2015 Difference" )
plt.xticks( [50,100,150,200,250] )
plt.yticks( [50,100,150,200,250] ) 

plt.tight_layout()
plt.savefig('/home/cparr/water_tracks/figs/drift_sets.png', dpi = 300)

In [34]:
img = src_dict_2012['actual snow depth']
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

In [35]:
rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()


/home/cparr/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:3: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  app.launch_new_instance()

In [36]:
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

In [37]:
rows, cols = img.shape
crow,ccol = rows/2 , cols/2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()


/home/cparr/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:6: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

In [ ]: