The Big Programming Thread - Page 777

#-- IMPORTS --#

import argparse                                                                         # This will import the argparse command-line parameters module.
import textwrap as _textwrap
import math                                                                             # This will import the math module
import numpy as np                                                                      # This will import the random number generation module
import matplotlib.pylab as plt

#-- FUNCTION DEFINITIONS --#

# Definition of the help text wrapping class:

class MultilineFormatter(argparse.HelpFormatter):
    def _fill_text(self, text, width, indent):
        text = self._whitespace_matcher.sub(' ', text).strip()
        paragraphs = text.split('|n ')
        multiline_text = ''
        for paragraph in paragraphs:
            formatted_paragraph = _textwrap.fill(paragraph, width, \
            initial_indent = indent, subsequent_indent = indent) + '\n\n'
            multiline_text = multiline_text + formatted_paragraph
        return multiline_text

#---------------------------------------------------------------------------------------        
def calculate_variance(x, y, x_foc, y_foc, \
                            x_min, x_max, y_min, y_max):

    # Calculation of the number of focal points:

    m = len(x_foc)                                                                      # Number of focal points = number of iterations in
                                                                                        # the first loop.
    # Calculation of the number of scales:

    scale = np.arange(1,45+1,1)                                                          # All scales.
    #scale = 1                                                                          # Only one scale.
    l = len(scale)

    # Preallocation of memory for all focal points:

    variance_avg_scale = np.zeros((m, 180))
    
    # Calculation of the number of analyzed points:
    
    n = len(x)
               
    # Calculation of transect properties:

    transect_angles = np.arange(0,180,1)                              
    slopes = np.tan((np.deg2rad(transect_angles)))
          
    # Selection of the focal point

    for j in range(m):
               
        # Creation of local coordinates /vectorized:

        i = np.arange(0,n,1)                                                                      # One iteration for each (x, y) point.
        x_local = x[i] - x_foc[j]
        y_local = y[i] - y_foc[j]

        # Limits of the studied area in local coordinates:
                
        x_min_local = x_min - x_foc[j]
        x_max_local = x_max - x_foc[j]
        y_min_local = y_min - y_foc[j]
        y_max_local = y_max - y_foc[j]
        
        # Translation to local polar coordinates:

        [angle, radius] = cart_to_pol(x_local,y_local)
        angle = np.rad2deg(angle)
        
        # Consolidation of coordinates data:

        data = np.vstack((angle, radius))

        # Preallocation of memory for all scales:

        transect_area = np.zeros(180)
        
        # Calculation of the transect area:
    
        for p in range(180):
        
            transect_angle = transect_angles[p]
            slope = slopes[p]
                
            transect_area[p] = calculate_transect_area(x_min_local, \
                                x_max_local, y_min_local, y_max_local, transect_angle, slope)
        
        # Preallocation of memory for all scales:

        variance_scale = np.zeros((l,180))

        # Selection of scale:

        for k in range(l):

            # Transect width:
            
            transect_width = scale
            
            # Preallocation of memory for all positions:

            wavelet_transform = np.zeros(180)
            variance_normalized = np.zeros(180)

            # Selection of angular position:

            for p in range(180):
               
                # Finding the points inside the specified transect:

                points_observed = observe_points(data, p, transect_width[k])

                # Extraction of the polar coordinates of the observed points:

                angle_observed = points_observed[0,:]
                radius_observed = points_observed[1,:]

                # Calculation of the number of observed points:

                observations = len(angle_observed)
                
                # Calculation of scaled wavelets for each point within the transect /vectorized:

                if observations == 0:
                    scaled_wavelet = [0]                                                # Has to be '[0]' instead of '0' because '0' cannot be used
                                                                                        # in sum (TypeError: 'int' object is not iterable).
                else:
                    o = np.arange(0,observations,1)
                    t = (angle_observed[o] - transect_angles[p])/scale[k]
                    scaled_wavelet = radius_observed*wavelet_function(t)
                
                wavelet_transform[p] = sum(scaled_wavelet)/scale[k]                     # No density measure.
                variance_normalized[p] = wavelet_transform[p]**2/transect_area[p]
            variance_scale[k,:] = variance_normalized
        variance_avg_scale[j,:] = sum(variance_scale)/l                                # All specified scales.
        #variance_avg_scale[j,:] = variance_scale                                       # Only one scale (1 degree).
    variance_avg_foc = sum(variance_avg_scale)/m
    return variance_avg_foc

#---------------------------------------------------------------------------------------
def cart_to_pol(x, y):

    angle = arctan(y, x)
    radius = np.sqrt(x**2 + y**2)
    return [angle, radius]

#---------------------------------------------------------------------------------------
def arctan(y, x):
    
    n = len(x);
    angle = np.zeros(n)
    for i in range(n):
        angle[i] = math.atan2(y[i], x[i])
        if angle[i] < 0:
            angle[i] = angle[i] + math.pi
    return angle

#---------------------------------------------------------------------------------------    
def observe_points(data, p, transect_width):

    # Finding the points inside the specified transect:

    log_ind = ( np.logical_or( np.logical_and( \
                               ((p - (transect_width/2)) <= data[0,:]), \
                               (data[0,:] <= (p + (transect_width/2))) ), \
                               np.logical_and( \
                               ((p - (transect_width/2) + 180) <= data[0,:]), \
                               (data[0,:] <= (p + (transect_width/2) + 180)) ) ) )
    
    points_observed = data[:,log_ind]
    return points_observed

#---------------------------------------------------------------------------------------    
def calculate_transect_area(x_min_local, x_max_local, y_min_local, y_max_local, \
                            transect_angle, slope):
                  
    # Slope variants:
    
    if (transect_angle == 0) or (transect_angle == 180):                                # Horizontal line.
        x_lim_1 = x_min_local
        x_lim_2 = x_max_local
        y_lim_1 = 0
        y_lim_2 = 0
    elif transect_angle == 90:                                                          # Vertical line.
        x_lim_1 = 0
        x_lim_2 = 0
        y_lim_1 = y_min_local
        y_lim_2 = y_max_local
    elif transect_angle > 0 and transect_angle < 90:
        
        # Boundary coordinates:

        if y_min_local/slope >= x_min_local:
            x_lim_1 = y_min_local/slope
            y_lim_1 = y_min_local
        else:
            x_lim_1 = x_min_local
            y_lim_1 = slope*x_min_local

        if y_max_local/slope <= x_max_local:
            x_lim_2 = y_max_local/slope
            y_lim_2 = y_max_local
        else:
            x_lim_2 = x_max_local
            y_lim_2 = slope*x_max_local

    else:
        
        # Boundary coordinates:

        if y_min_local/slope >= x_max_local:
            x_lim_1 = x_max_local
            y_lim_1 = slope*x_max_local
        else:
            x_lim_1 = y_min_local/slope
            y_lim_1 = y_min_local

        if y_max_local/slope <= x_min_local:
            x_lim_2 = x_min_local
            y_lim_2 = slope*x_min_local
        else:
            x_lim_2 = y_max_local/slope
            y_lim_2 = y_max_local
    
    transect_area = (x_lim_1)**2 + (y_lim_1)**2 + (x_lim_2)**2 + (y_lim_2)**2
    return transect_area
                                                                                        # This is not real area of the transect.
                                                                                        # It's area divided by omega/2, but that holds
                                                                                        # for all transects and omega is constant.    

#---------------------------------------------------------------------------------------
def wavelet_function(t):

    mexican_hat = 2/(3**(1/2))*math.pi**(-1/4)*(1 - 4*t**2)*np.exp(-2*t**2)            # Wavelet function normalized to have unit energy.
    return mexican_hat

#-- PARAMETERS --#

parser = argparse.ArgumentParser(description = 'This program quantifies the degree \
of banding in the microstructure by means of angular wavelet analysis. The only \
valid name of the input data is "coordinates.txt". \n',
formatter_class = MultilineFormatter)

# Optional arguments:

parser.add_argument("-p", "--edge_parameter",                                           # Edge parameter limiting the focal points area.
                    type = float,
                    default = 0.25,
                    help = 'Edge parameter determines the boundaries of the focal \
points area, as measured from the edge of the analyzed image. The default value of \
the edge parameter is 0.25.')

parser.add_argument("-f", "--focal_points",                                             # Shows points within the focal points area.
                    action = 'store_true',
                    help = 'Shows points within the focal points area.')

args = parser.parse_args()

#-- IMPORT OF DATA --#

D = np.loadtxt("coordinates.txt", delimiter = ' ')

#-- EXTRACTION OF COORDINATES --#

x_extracted = D[:,0]
y_extracted = D[:,1]

# Transposition into vectors:

x = x_extracted.T
y = y_extracted.T

# Verification of vector lengths:

if len(x) != len(y):
    print('The imported data is corrupted: the numbers of the x and y coordinates \
do not match')

#-- BOUNDARY CONDITIONS --#

# Limits of the analyzed space:

x_min = min(x)
x_max = max(x)
y_min = min(y)
y_max = max(y)

# Edge region parameters for each coordinate:

edge_parameter = args.edge_parameter

edge_parameter_x = edge_parameter*(x_max - x_min)
edge_parameter_y = edge_parameter*(y_max - y_min)

# Limits of the focal point space:

x_foc_min = x_min + edge_parameter_x
x_foc_max = x_max - edge_parameter_x
y_foc_min = y_min + edge_parameter_y
y_foc_max = y_max - edge_parameter_y

#-- EXTRACTION OF FOCAL POINTS --#

# Generation of focal points:

x_foc = x[np.logical_and(np.logical_and(x > x_foc_min,x < x_foc_max), \
        np.logical_and(y > y_foc_min,y < y_foc_max))]
y_foc = y[np.logical_and(np.logical_and(x > x_foc_min,x < x_foc_max), \
        np.logical_and(y > y_foc_min,y < y_foc_max))]

# Focal points plot:

if args.focal_points:

    fig = plt.figure()
    ax1 = fig.add_subplot(111)

    ax1.scatter(x, y, s = 1, color = 'black')                                           # All points.
    ax1.scatter(x_foc, y_foc, s = 1, color = 'red')                                     # Focal points.
    plt.show()

#-- ANALYSIS --#

overall_variance = calculate_variance(x, y, x_foc, y_foc, x_min, x_max, y_min, y_max)

#-- OUTPUT --#

# Original variance peak:

print('Original variance peak:')
print(max(overall_variance))

# Plot:

plt.plot(np.arange(0,180,1), overall_variance)
plt.show()