%matplotlib inline import numpy as np from skimage import io import matplotlib.pyplot as plt class stats: def __init__( self, head = None, tail = None, url = None, normalize = True, periodic = np.array([]), shift = True, display = True, cutoff = None): self.periodic = np.array( periodic ) if url is None: self.head = head self.tail = tail elif url is not None: # URL is a quickstart option to compute statistics really fast out of the box self.head = LoadImage( url ) self.tail = None self.url = url # Initialize cutoff if cutoff is None: cutoff = SetCutoff( self.head.shape ) elif isinstance(cutoff, int): cutoff = np.array( [cutoff] ) elif isinstance(cutoff, list): cutoff = np.array( cutoff ) self.cutoff = cutoff self.Stats( normalize, shift, display, cutoff ) def Stats( self, normalize, shift, display, cutoff ): sz = Periodic2Size( self.head.shape, self.periodic) self.spatial = convolve( self.head, self.tail, sz ) vectors = AppendIndex( self.head.shape, self.spatial ) self.vectors = vectors # self.units.values defines the normalized vector values if normalize: # ``normalize`` passed from __init__ method, denom = self.__Normalize() self.denom = denom self.spatial = np.divide( self.spatial, denom ) # slice/truncate statistics for ax in np.arange( len( self.spatial.shape ) ): # This indexing conditions on cutoff could be problematic. Need size checking cutoffid = min( cutoff.size - 1, ax) # subtract 1 because ``ax`` starts at zero self.spatial = np.compress( np.abs( self.vectors[ax] ) < cutoff[ cutoffid ], self.spatial, axis = ax ) self.vectors[ax] = np.compress( np.abs( self.vectors[ax] ) < cutoff[ cutoffid ] , self.vectors[ax] ) if shift: # FFT Shift self.spatial = np.fft.fftshift( self.spatial ) # Shift indices for ii, ax in enumerate( self.vectors ): self.vectors[ii] = np.fft.fftshift( ax ) # Switch X,Y # http://stackoverflow.com/questions/24791614/numpy-pcolormesh-typeerror-dimensions-of-c-are-incompatible-with-x-and-or-y plt.pcolor(self.vectors[1], self.vectors[0], self.spatial.real) # Requires special treatment after slicing def numer( self ): # Determine denominator vector axes vectors = AppendIndex( self.head.shape, self.denom ) # trim the denominator denom = self.denom for ax in np.arange( len( denom.shape ) ): # This indexing conditions on cutoff could be problematic. Need size checking cutoffid = min( self.cutoff.size - 1, ax) # subtract 1 because ``ax`` starts at zero denom = np.compress( np.abs( vectors[ax] ) < self.cutoff[ cutoffid ], denom, axis = ax ) # Unshifted matrices will have vectors whose first index is zero if not np.all( [v[0]==0 for v in s.vectors] ): denom = np.fft.fftshift( denom ) return np.multiply( self.spatial, denom ) def __Normalize( self ): # Decision on the denominator here sz = Periodic2Size( self.head.shape, axes = self.periodic) if len( self.periodic ) >= 0: # Swap this to a leaner function method = 'convolve' denom = convolve( head = np.ones( shape = self.head.shape ), tail = None, sz = sz) else: print len( self.periodic ) return method, denom def npNormalization( sz, dims ): # dims is a d x 2 array # the first and second columns are the upper and lower bounds of the vector # It computes something requiring the manhattan return def LoadImage( url ): A = np.array( io.imread( url ) ) # Process the Image A = Binarize( A ) return A def convolve( head, tail = None, sz = None ): # Name the variables based on the reference and moving features if sz is None: sz = np.array( head.shape ) fA = np.fft.fftn( head, s = sz ) if tail is None: # Use [Wiener Kichin](http://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem) fA = np.abs( fA ) fA = fA ** 2 else: # Brute Force FFT approach for different states fA = np.multiply( fA , np.conj( np.fft.fftn( tail, s = sz ) ) ) fA = np.fft.ifftn( fA, s = sz ) fA = np.real_if_close( fA ) return fA def AppendIndex( osz, fA ): # osz - Original Size of the raw data - A.shape # Define the unit vector distances for the statistics nsz = fA.shape ind = list() for id, ax in enumerate( nsz ): a = np.arange(0, ax) a[ a > osz[id]/2 ] = a[ a > osz[id]/2 ] - ax ind.append( a ) # A list of unit vector indices for each axis return ind def SetCutoff( sz ): # Set the cutoff criteria for the vectors return np.floor(np.array( sz )/2) def Periodic2Size( oshape, axes, mult = .6 ): # Determine the size of the FFT domain for the convolution # Initialize the size of nonperiodic Statistics FFT size nshape = np.array( oshape ) + np.array( oshape ) * mult if len( axes ) > 0: # Reassign periodic axes where they exist for ax in axes: nshape[ax] = oshape[ax] return nshape.round().astype('int') def Binarize( A ): # Binarize a grayscale image downloaded from the web return np.round( A.astype( 'double' ) /255 ) A = np.zeros( (5,5)) A = np.eye( 5, 5 ) s = stats( head = A ) print s.numer() s = stats( head = A ) s.numer().real url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' A = io.imread( url ) A = Binarize( A ) period = [0,1] s = stats( head = A, periodic=period, display=False ) print "Volume Fractions:: %f / %f" % (s.spatial.real.max(), s.head.mean()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) # plt.imshow( s.numer().real ) # print "Counts:: %f / %f" % (s.numer().real.max(), s.head.sum()) # url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' # A = io.imread( url ) period = [0,1] s = stats( head = A, shift = True, normalize=True, periodic=period ) print "Volume Fractions:: %f / %f" % (s.spatial.real.max(), s.head.mean()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) #url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' #A = io.imread( url ) #A = Binarize( A ) B = 1-A period = [0,1] s = stats( head = A, tail=B, periodic=period ) print "Minimum Magnitude:: %f" % (np.abs(s.spatial.real).min()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) # url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' # A = io.imread( url ) B = 1-A period = [0,1] s = stats( head = A, tail=B, periodic=period, cutoff = [100, 50] ) print "Minimum Magnitude:: %f" % (np.abs(s.spatial.real).min()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) B = 1-A period = [0,1] s = stats( head = A, tail=B, periodic=period, cutoff = 50 ) print "Minimum Magnitude:: %f" % (np.abs(s.spatial.real).min()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) #url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' #A = io.imread( url ) #A = Binarize( A ) B = 1-A s = stats( head = A, tail=B ) print "Minimum Magnitude:: %f" % (np.abs(s.spatial.real).min()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) #url = 'https://farm8.staticflickr.com/7368/12972465885_1bcb9717bd_c.jpg' #A = io.imread( url ) #A = Binarize( A ) s = stats( head = A ) print "Volume Fractions:: %f / %f" % (s.spatial.real.max(), s.head.mean()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) period = [1] s = stats( head = A, periodic=period ) print "Volume Fractions:: %f / %f" % (s.spatial.real.max(), s.head.mean()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) s = stats( url = url ) print "Volume Fractions:: %f / %f" % (s.spatial.real.max(), s.head.mean()) print "Maximum Imaginary Value :: " + str( s.spatial.imag.max() ) # Parametric Analysis of Imaginary Components After Communiting Convolution Functions With np.fft for v in np.arange( 101,1002,100): sz = np.array( (v,v) ) A = np.ones( shape=sz ) B = np.concatenate(( np.concatenate( (A, np.zeros( sz )), axis = 0 ), np.concatenate( (np.zeros( sz ), np.zeros( sz )), axis = 0 ) ),axis = 1) fB = np.fft.fft2(B) fB = fB * np.conj( fB ) fB = np.fft.ifft2( fB ) fA = np.fft.fft2( A, s = 2 * sz ) fA = fA * np.conj( fA ) fA = np.fft.ifft2( fA ) print v,fA.imag.max(),fB.imag.max()