The Big Programming Thread - Page 768

clc
clear

%%% PARAMETERS %%%

answer = input('Should the default value be used for the edge region parameter? Y/N [Y]: ','s');
if isempty(answer)                                                         % If the user simply presses 'Enter', the answer is assumed to be 'Y'.
    answer = 'Y';
end

disp(' ')                                                                  % Empty line for clarity.

if strcmp(answer, 'Y') == 1
    edge_parameter = 0.25;
    disp('The default value of the edge region parameter (0.25) will be used.')
else
    edge_parameter = input('Please, specify the value of the edge region parameter: ');
end

%%% IMPORT OF DATA %%%

filename = 'banding_coordinates.txt';
delimiterIn = ' ';
headerlinesIn = 1;
D = importdata(filename,delimiterIn,headerlinesIn);

%%% EXTRACTION OF COORDINATES %%%

x_extracted = D.data(:,1);
y_extracted = D.data(:,2);

% Transposition into vectors:

x = transpose(x_extracted);
y = transpose(y_extracted);

% Verification of vector lengths:

if length(x) ~= length(y)
    disp('The imported data is corrupted: the numbers of the x and y coordinates do not match')
end

%%% 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_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:

[in,on] = inpolygon(x, y, [x_foc_min, x_foc_max], [y_foc_min, y_foc_max]);

x_foc = x(in);
y_foc = y(in);

% Focal points plot:

answer = input('Would you like to display the generated focal points? Y/N [N]: ','s');
if isempty(answer)                                                         % If the user simply presses 'Enter', the answer is assumed to be 'Y'.
    answer = 'N';
end

disp(' ')                                                                  % Empty line for clarity.

if strcmp(answer, 'Y') == 1
    figure

    scatter(x,y,'.')                                                       % all points
    axis equal

    hold on
    scatter(x_foc,y_foc,'r+')                                              % focal points
    hold off
end

%%% ANALYSIS %%%
tic
overall_variance = calculate_variance_test6_3(x, y, x_foc, y_foc, ...
                   x_min, x_max, y_min, y_max);
toc
%%% OUTPUT %%%

% Original variance peak:

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

% Plot:

figure

plot(1:1:180, overall_variance)

function variance_avg_foc = calculate_variance_test6_3(x, y, x_foc, y_foc, ...
                            x_min, x_max, y_min, y_max)

    % Calculation of the number of focal points:

    m = length(x_foc);                                                     % Number of focal points = number of iterations in the first loop.

    % Calculation of the number of scales:

    scale = 1:1:45;                                                       % All scales.
    l = length(scale);

    % Preallocation of memory for all focal points:

    variance_avg_scale = zeros(m, 180);
    
    % Calculation of the number of analyzed points:
    
    n = length(x) - 1;
                
    % Calculation of transect properties:

    p = 1:1:180;                              
    transect_angles = p - 1;
    slopes = tand((transect_angles));
          
    % Selection of the focal point

    for j = 1:m
        
        % Removal of the focal point from the pool of analyzed points:    
        
        index = find( (x(1,:) == x_foc(j)) & (y(1,:) == y_foc(j)) );
                
        if length(index) ~= 1                                              % If there is more than once point with the x_foc(j), y_foc(j)
            index = index(1);                                              % coordinates, only one point will be subtracted from the future
        end                                                                % local coordinates pool.
        
        x_without_foc = x(1:end ~= index);
        y_without_foc = y(1:end ~= index);
        
        % Creation of local coordinates /vectorized:

        i = 1:1:n;                                                        % One iteration for each (x, y) point.
        x_local = x_without_foc(i) - x_foc(j);
        y_local = y_without_foc(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 = rad2deg(angle);
        
        % Consolidation of coordinates data:

        data = [angle; radius];

        % Preallocation of memory for all scales:

        transect_area = zeros(1, 180);
        
        % Calculation of the transect area:
    
        for p = 1: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);
        end
        
        % Preallocation of memory for all scales:

        variance_scale = zeros(l, 180);

        % Selection of scale:

        for k = 1:l

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

            wavelet_transform = zeros(1, 180);
            variance_normalized = zeros(1, 180);

            % Selection of angular position:

            for p = 1: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(1,:);
                radius_observed = points_observed(2,:);

                % Calculation of the number of observed points:

                observations = length(angle_observed);
                
                % Calculation of scaled wavelets for each point within the transect /vectorized:

                if observations == 0
                    scaled_wavelet = 0;
                else
                    o = 1:1:observations;
                    t = (angle_observed(o) - transect_angles(p))./scale(k);
                    scaled_wavelet = radius_observed.*wavelet_function(t);
                end
                
                wavelet_transform(p) = sum(scaled_wavelet)/scale(k);
                variance_normalized(p) = wavelet_transform(p)^2/transect_area(p);
            end
            variance_scale(k,:) = variance_normalized;
        end
        variance_avg_scale(j,:) = sum(variance_scale)./l;                 % All specified scales.
    end
    variance_avg_foc = sum(variance_avg_scale)./m;
end

%--------------------------------------------------------------------------
function [angle, radius] = cart_to_pol(x, y)

    angle = arctan(y, x);
    radius = hypot(x,y);

end

%--------------------------------------------------------------------------
function angle = arctan(y, x)
    
    n = length(x);
    angle = zeros(1, n);
    parfor i = 1:n                                                         % /vectorized
        angle(i) = atan2(y(i), x(i));
        if angle(i) < 0
            angle(i) = angle(i) + pi;
        end
    end
    
end

%--------------------------------------------------------------------------
function points_observed = observe_points(data, p, transect_width)

    % Finding the points inside the specified transect:

    step_size = 1;

    log_ind = ( (p - 1 - (transect_width/2))*step_size <= data(1,:) & ...
                 data(1,:) <= (p - 1 + (transect_width/2))*step_size ) | ...
              ( ((p - 1 - (transect_width/2))*step_size + 180) <= data(1,:) & ...
                 data(1,:) <= ((p - 1 + (transect_width/2))*step_size + 180) );
    points_observed = data(:,log_ind);

end

%--------------------------------------------------------------------------
function transect_area = calculate_transect_area(x_min_local, ...
                         x_max_local, y_min_local, y_max_local, transect_angle, slope)
                  
    % Slope variants:
    
    if (transect_angle == 0) || (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;
    elseif 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;
    elseif transect_angle > 0 && 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;
        end

        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;
        end
    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;
        end

        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;
        end
    end
    
    transect_area = (x_lim_1)^2 + (y_lim_1)^2 + (x_lim_2)^2 + (y_lim_2)^2; % 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.
end

%--------------------------------------------------------------------------
function mexican_hat = wavelet_function(t)                                 % /vectorized

    mexican_hat = 2./(sqrt(3)).*pi.^(-1/4).*(1 - 4.*t.^2).*exp(-2.*t.^2);  % Wavelet function normalized to have unit energy.
                                                                           % Vectorizing with "vectorize()" does not work.
                                                                           % It cannot tell if I want to use ^ for matrices (mpower)
                                                                           % or scalars (power).
    
end

points_observed = observe_points(data, p, transect_width(k));

t = (angle_observed(o) - transect_angles(p))./scale(k);
scaled_wavelet = radius_observed.*wavelet_function(t);

% Calculation of scaled wavelets for each point within the transect /vectorized:

if observations == 0
    scaled_wavelet = 0;
else
    log_ind = (radius_observed ~= 0);
    %o = 1:1:observations;
    scaled_wavelet = radius_observed(log_ind).*wavelet_function((angle_observed(log_ind) - transect_angles(p))./scale(k));
end