% Authors: Steven Su and Robert Kim
% Modified from Existing Code from Nathanael Kuo
% Date: May 17, 2011
% Interventional Photoacoustic Registration
% 600.446: CIS II, Spring 2011
% The Johns Hopkins University

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gets the centroids of the photoacoustic signal in the ultrasound images.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [us_cx us_cy Bpost] = getUltrasoundPoint(path)
close all;

nCh = 128;              % number of channels
reRoute = true;         % true: transducer element (correct image), false: DAQ element
chanls = ones(1,nCh);   % what channels to read (DAQ element), for each channel set to 1
fulAprSz = 64;          % aperature size

% Sets correct directory
D = dir(path);
if nargin == 1
    idx = find([D.isdir]);
    folder = D(idx(end)).name;
end

% Initializes x,y coordinates for the centroid of the US image
us_cx = 0;
us_cy = 0;

% Loops through all 5 frames of each US image. Calculates the weighted
% centroid of each frame for the PA signal. Peforms some image processing
% algorthms to threshold image. 
for frameN = 1:5
    [hdr, RF] = readDAQ(path, chanls, frameN, reRoute); % reads DAQ data
    Bpre = sqrt( abs( hilbert(RF) ) ) ; % Computes prebeamform image
    
    out = beamformDAQ(RF, fulAprSz); % Computes beamformed image
    Bpost = sqrt( abs( hilbert(out) ) ) ; % Processes beamformed image
    
    Bpost2=Bpost(600:end,:); % Temporarily removes noise bar that occurs in images for image processing purposes
    mean_Bpost = mean(Bpost2,2);% Calculates the mean intensity of the image 
    third_percentile = prctile(mean_Bpost,99.9); % Calculates the 99.9 percentile intensity of the image, user can change this variable.
    
    % Thresholds image using the selected percentile above.
    threshold_image = Bpost2;
    for i=1:size(Bpost2,1)
        for j=1:size(Bpost2,2)
            if(Bpost2(i,j)<third_percentile)
                threshold_image(i,j) = 0;
            else
                threshold_image(i,j) = 1;
            end
        end
    end
    
    % Another thresholding algorthm. Removes any segmented regions with
    % less then num_pix pixels.
    if(strfind(path,'us_data_low_right_corner'))
        num_pix = 300;
    else
        num_pix = 500;
    end
    bw = threshold_image;
    bw2 = bwareaopen(bw, num_pix,18);
    
    % Replaces the previously removed band of pixels where the noise
    % artifact was in orginal image with black pixels.
    bw3=bw2;
    Bpost4=[zeros(600,size(bw3,2)); bw3];
    
    % Calculates the centroid of the thresholded image
    L = bwlabel(Bpost4);
    s = regionprops(L, 'Centroid');
    v2=cat(1,s.Centroid);
    cx2=v2(end,1);
    cy2=v2(end,2);
    
    % Sums the cumulative x and y centroid coordinates for computing the
    % average later on.
    us_cx = us_cx + cx2;
    us_cy = us_cy + cy2;
    
    pause(2);      

    close all;    
end

% Computes the mean centroid for the image based on all the frames.
us_cx = us_cx/5;
us_cy = us_cy/5;

