

% function R = read_surface_extract(day,REstar,Bo,time)
%
%  reads in the raw surface extract file for a particular simulation and time

function [Xs,Ys,Rs] = read_processed_surface_extract(day,REstar,Bo,time)


% this function below reads teh raw surface extracts, no filtering and no sorting

R = read_surface_extract(day,REstar,Bo,time);

RR = R(:,3:4);
[A,igood] = unique(RR,'rows','stable');
R = R(igood,:);
nR = size(R,1);
nRR = size(RR,1);



[nR,mR] = size(R);
PC = [R(:,3:4)  0*R(:,3)];
RAWpc = pointCloud( PC );

Xraw = R(:,3);  Yraw = R(:,4);


denoise_fac = 4;
[Rpc, good_index] = pcdenoise(RAWpc,'NumNeighbors',2,'Threshold',denoise_fac);
numRpc = size(Rpc.Location,1);
Rdn = R(good_index,:);   % same as Rpc but now keeps only denoised points

      [labels,nclusters]=pcsegdist(Rpc,0.01);
      ncpoints = zeros(nclusters,1);
     ii = [];
     for k=1:nclusters,
         ncpoints(k) = sum( labels==k);
         if ncpoints(k)>1,
             itmp = find(labels==k);
             ii = [ii; itmp];
         end;
     end;

     Rsd = Rdn(ii,:);  % now segmented distributed
     nRpc = select(Rpc,ii);
     nnRpc = size(nRpc.Location,1);
     [Xn,Yn,ds,Rs]=sort_surface_extract(nRpc.Location(:,1), nRpc.Location(:,2),Rsd);
     nxn = length(Xn);   C=[1:nxn]/nxn';
     %ifar = find(Xn>46);
     %if (isempty(ifar))
     %    disp('bad: ifar is empty')
     %    ifar = find(Xn>43);
     %end;
     %min_ifar = min(ifar);
     Xs = Xn;   Ys = Yn;
     ns = size(Xs,1);
     disp(sprintf('Restar=%d Bo=%d,  time=%.2f, # Raw points = %d,  # Unique Points = %d',REstar,Bo,time,nRR,nR))
     disp(sprintf('denoise_fac=%f;  denoise numR=%d, nclusters=%d, pcsegdist numR=%d; ns=%d',denoise_fac,numRpc,nclusters,nnRpc,ns))     

     %     fprintf (fp_surface_file, '%g %g %g %g %g %g %g %g\n', x, y, xpos[], ypos[], u.x[], u.y[], vort[], p[]);

%      figure(1)
%      plot(Xraw,Yraw,'k.',Rpc.Location(:,1),Rpc.Location(:,2),'bo')
%      hold on
%     plot(Xs,Ys,'rx')
%      keyboard
          return
      
    
