function so = average_p(si, whichway,bounds,npts)
%SAXS/AVERAGE_P Average SAXS image intensities in pixel-by-pixel mode
%   AVERAGE(S, 'x',points) returns a SAXS image object S with image
%   columns collapsed and replaced by the average intensity in that column.
%   It uses the original 2D data and bins all the pixels into their
%   respective positions to create a 1D data set.
%
%   In other words, AVERAGE_raw returns the one-dimensional curve of average
%   image intensity versus whatever the x coordinate is, usually momentum
%   transfer (image radius).
%
%   This routine differs from average_raw in that the reductions are
%   here performed on a pixel-by-pixel basis.

if isempty(si.reduction_state) %defining reduction state when no reductions
    si.reduction_state='000000000000';
end

%setting up some defaults
if nargin == 1
    whichway = 'x';  % default
    npts=200; %default
elseif nargin == 2
    npts=200; %default
elseif ~ (isequal(whichway, 'x') || isequal(whichway, 'y'))
    error('The second argument must be ''x'' or ''y''.')
end

error_norm_fact=si.error_norm_fact;



warning off MATLAB:divideByZero
si=insert_saxs_pixels(si,si.raw_pixels); %inserting raw_pixels to be reduced
so=pix2pix_red(si); % that is the reduced data
%This hack is included to fix a bug when the user has not entered any calibrations.
if isempty(si.pixelcal)
    si.kcal=1.;
    si.pixelcal=[1.0 1.0];
    si.koffset=0;
    si.wavelength=[];
else
    %this is a correction to allow for WAXS si.kcal becomes
    %deltaq_at_0*si.pixelcal
    si.kcal=tan(2*asin(si.kcal*si.wavelength/4/pi))/2*4*pi/si.wavelength;
end
% now for the averaging
switch whichway
    case 'x'
        [x,yave,yrelerr,yrelerr2]=azimbinning(so.pixels,bounds(1),bounds(1)+bounds(3),npts,...
            so.xaxisbintype,bounds(2),bounds(2)+bounds(4),so.raw_center(2),so.raw_center(1),...
            si.kcal/si.pixelcal(2),si.kcal/si.pixelcal(1),si.wavelength,error_norm_fact);
        % now make an intensity correction to bring the units into
        % intensity per solid angle (mr^2)
        if isempty(si.wavelength)
            solid_angle=1;
        else
            solid_angle=(si.kcal^2/si.pixelcal(1)/si.pixelcal(2))*(si.wavelength/2/pi)^2*...
                cos(2*asin(x'*si.wavelength/4/pi)).^3*1E6;
        end
        
        yave=yave./solid_angle;
        switch so.type
            case 'xy'
                so.type = 'x-ave';
                event = 'average vs x';
            case 'rth'
                so.type = 'r-ave';
                event = 'average vs r';
        end
        so.pixels=yave';
        so.relerr=yrelerr';
        so.xaxisfull=x;
    case 'y'
        % we use a different approach here due to problems with discretization
        % issues at odd angles
        % here we use the pol2cart method of the original SAXSGUI versions.
        %
        % also reduction is not yet implemented here
        % first we transform to polar


        st=calc_azim_int(so,bounds,npts); %calculating the average vs azimuthal angle

        %this following calculation is valid for all following
        %averages...but needs to be done at least once.
        len = length(st.pixels);
        x = st.xaxis;
        y = st.yaxis;
        x = [x(1) : (x(2) - x(1)) / (len - 1) : x(2)];
        y = [y(1) : (y(2) - y(1)) / (len - 1) : y(2)];

        % calibrated
        if si.kcal
            UnitsPerPixel=stemp.kcal/mean(st.pixelcal);
        else
            UnitsPerPixel=1;
        end

        dx=x(2)-x(1);
        dy=y(2)-y(1);
        xmin=x(1);
        ymin=y(1);
        xmax=x(len);
        ymax=y(len);

        NumOfPix=pi*(xmax^2-xmin^2)/UnitsPerPixel^2*(dy)/360;

        % now we have what we need for calculating the average and errors
        % on the actual data.
        yave=st.pixels;
        yrelerr=1./sqrt(NumOfPix.*st.pixels)*error_norm_fact;

        so.type = 'th-ave';
        event = 'average vs theta';
        so.pixels=yave;
        so.relerr=yrelerr;
        so.xaxisfull=y;
end

warning on MATLAB:divideByZero
so.history = strvcat(so.history, event);
so.xaxis=[bounds(1) bounds(1)+bounds(3)];
so.yaxis=[bounds(2) bounds(2)+bounds(4)];




