function [fos]=first_order_statistics(volume_vec)
%计算三维体的一阶统计特征
%（记得在使用前检测最小值是否为0，DICOM数据截断是否为-1024）
%（过程中直方图的首个值为最低像素值的计数。
%  因为即使排除了0值，空气中仍有小值对结果产生影响；
%  而对于肿瘤分割和roi框定的数据，其值应为所涉及区域内的真实值)
%输入为体数据的数据向量

fos = struct('energy',{0.0}, 'entropy',{0.0}, 'krutosis',{0.0},...
    'maximum', {0.0},'mean',{0.0},'mean_absulute_deviation', {0.0},'median',...
    {0.0}, 'minimum',{0.0}, 'range',{0.0}, 'root_mean_square',{0.0},'skewness',...
    {0.0}, 'standard_deviation', {0.0},'uniformity', {0.0},'variance',{0.0},...
    'mass',{0.0},'entropy_p',{0.0},'uniformity_p',{0.0});

%% 生成数据的灰度直方图
volume_max=max(volume_vec);
volume_min=min(volume_vec);
volume_hist0=hist(volume_vec,(volume_max-volume_min+1));
% volume_hist=zeros(volume_max+1,1);
% volume_hist((volume_min+1):(volume_max+1))=volume_hist0;
n=size(volume_vec,1);
volume_hist=volume_hist0;
volume_hist_p=volume_hist./n;

% if volume_min>=0
%     t=0;
% else
%     t=1;
% end
% for i=1:1:z
%     v=volume(:,:,i);
%     i_max=max(v(:));
%     i_min=min(v(:));
%     if i_max<=0
%         continue
%     end
%     h=hist(v,(i_max-i_min+1));
%     h_t=sum(h');
%     h=h_t';
%     if i_min==0
%         volume_hist(1:i_max)=volume_hist(1:i_max)+h(2:1:end);
%     elseif i_min>0
%         volume_hist(i_min:i_max)=volume_hist(i_min:i_max)+h;
%     elseif i_min<0
%         volume_hist(1:i_max)=volume_hist(1:i_max)+h((2-i_min):1:end);
%     end
% end


%% 提取特征
log_resdue=0.000001;
volume_pixelnum=size(volume_vec,1);
volume_hist_log2=log2(volume_hist+log_resdue);
volume_hist_p_log2=log2(volume_hist_p+log_resdue);

fos.mean = sum(volume_vec)/volume_pixelnum;
fos.energy = sum(volume_vec.^2);
fos.entropy = sum(volume_hist.*volume_hist_log2);
fos.krutosis = (sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^4)/volume_pixelnum)...
    /(sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^2)/volume_pixelnum).^2;
fos.maximum = max(volume_vec);
fos.mean_absulute_deviation =sum(abs(volume_vec-ones(volume_pixelnum,1).*fos.mean))/volume_pixelnum;
fos.median = median(volume_vec);
fos.minimum = min(volume_vec);
fos.range = fos.maximum-fos.minimum;
fos.root_mean_square = sqrt(sum(volume_vec.^2)/volume_pixelnum);
fos.skewness = (sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^3)/volume_pixelnum)...
    /power(sqrt(sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^2)/volume_pixelnum),3);
fos.standard_deviation = sqrt(sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^2)/(volume_pixelnum-1));
fos.uniformity = sum(volume_hist.^2);
fos.variance = sum((volume_vec-ones(volume_pixelnum,1).*fos.mean).^2)/(volume_pixelnum-1);
fos.mass=sum(volume_vec);
fos.entropy_p=sum(volume_hist_p.*volume_hist_p_log2);
fos.uniformity_p=sum(volume_hist_p.^2);
