function [out]=peakInt(inmat,verpts,horpts,wnd,refn)
% *** intensity at 9 pts around peak-top ***
% Y. Matsuki, August, 2010
% usage:
% [out]=ppm2pts(in,sw,np,cent,freq,dir)
% =====================================

npk=length(verpts(:,1));
verptsmod=verpts;
horptsmod=horpts;

int=zeros(npk,1);
nr=length(inmat(:,1));
nc=length(inmat(1,:));

vwidth=1; % arbitrary!
% hwidth=4; % arbitrary!
% hwidth=2; % arbitrary!
hwidth=1; % arbitrary!
for m=1:npk
    if refn>0 % added, aug1310.
        if verpts(m)-vwidth<1
            a=1;
            b=verpts(m)+vwidth;
        elseif verpts(m)+vwidth>nr;
            b=nr;
            a=verpts(m)-vwidth;
        else
            a=verpts(m)-vwidth;
            b=verpts(m)+vwidth;
        end
        
        if horpts(m)-hwidth<1
            c=1;
            %         disp('1')
            d=horpts(m)+hwidth;
        elseif horpts(m)+hwidth>nc;
            d=nc;
            %         disp('2')
            c=horpts(m)-hwidth;
        else
            c=horpts(m)-hwidth;
            d=horpts(m)+hwidth;
        end
        
        temp=real(inmat(a:b,c:d));
        [c,i]=max(temp);
        [cc,ii]=max(c);
        %     temp(i(ii),ii); % max value.
        
        vsize=length(temp(:,1));
        vshift=ceil((vsize-i(ii))-(vsize/2));
        verptsmod(m)=verpts(m)-vshift;
        
        %     if i(ii)<4
        %         verptsmod(m)=verpts(m)-(4-i(ii));
        %     elseif i(ii)>4
        %         verptsmod(m)=verpts(m)+(i(ii)-4);
        %     end
        
        hsize=length(temp(1,:));
        hshift=ceil((hsize-ii)-(hsize/2));
        horptsmod(m)=horpts(m)-hshift;
        
        %     if ii<4
        %         horptsmod(m)=horpts(m)-(4-ii);
        %     elseif ii>4
        %         horptsmod(m)=horpts(m)+(ii-4);
        %     end
    end
    
    
    %% collect intensities.
    % wnd=1;
    
    % trim intensity summation window.
    a=verptsmod(m)-wnd;
    b=verptsmod(m)+wnd;
    c=horptsmod(m)-wnd;
    d=horptsmod(m)+wnd;
    int(m)=0;
    % vertical axis
    if a < nr/2 &&  b > nr/2
        pts=b-nr/2; % sigf2(m,1) is 0 or negative.
        b=nr/2;
        aa=nr-pts+1;
        bb=nr;
        cc=c;
        dd=d;
        int(m)=int(m)+sum(sum(real(inmat(aa:bb,cc:dd))));
    elseif b > nr
        pts=b-nr; % sigf2(m,2) is larger than nc.
        b=nr;
        aa=nr/2-pts+1;
        bb=nr/2;
        cc=c;
        dd=d;
        int(m)=int(m)+sum(sum(real(inmat(aa:bb,cc:dd))));
    end
    
    % horizontal axis
    if c < 1
        pts=1-c; % c is 0 or negative.
        c=1;
        cc=nc-pts+1;
        dd=nc;
        aa=a;
        bb=b;
        int(m)=int(m)+sum(sum(real(inmat(aa:bb,cc:dd))));
    elseif d > nc
        pts=d-nc; % pts>=1.
        d=nc;
        cc=1;
        dd=pts;
        aa=a;
        bb=b;
        int(m)=int(m)+sum(sum(real(inmat(aa:bb,cc:dd))));
    end
    
    int(m)=int(m) + sum(sum(real(inmat(a:b,c:d))));
    %         int(m)=sum(sum(real(inmat(verptsmod(m)-wnd:verptsmod(m)+wnd,horptsmod(m)-wnd:horptsmod(m)+wnd))));
end

out=[verptsmod horptsmod int];











