function Cat=sim_seg_moments(Sim,SimSeg,varargin)
%--------------------------------------------------------------------------
% sim_seg_moments function                                          ImPhot
% Description: Compute the segments' moments for a set of structure images.
% This function requires a pre-knowledge of the segmentation image, which
% may be computed by the sim_segmentation function.
% Input  : - Set of images (See images2sim.m for options).
%          - Set of segmentation images (See images2sim.m for options).
%            The number of segmentation images must be equal to the number
%            of images or to one.
%            A segmentation image contains the value 1 in the background,
%            and the integer value N for the N-th segment (e.g., source).
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'OutPar' - A cell containing the parameters that will be
%                       computed for each segment.
%                       Default is {'Iseg','X','Y','X2','Y2','XY','ALPHA_J2000','DELTA_J2000','NPIX',...
%                       'X_PEAK','Y_PEAK','ALPHA_PEAK','DELTA_PEAK','FLUX_PEAK','FLUX'}
%                       Where:
%                       'X','Y' - first moment in image coordinates.
%                       'X2','Y2','XY' - second moments in image coordinates.
%                       'ALPHA_J2000','DELTA_J2000' - first moment in 
%                                 J2000.0 RA/Dec.
%                       'NPIX' - Number of pixels in segment.
%                       'X_PEAK','Y_PEAK' - Position of the pixel with
%                                 the maximum intensity.
%                       'ALPHA_PEAK','DELTA_PEAK' - Position of the pixel
%                                 with maximum intensity, in celestial
%                                 coordinates.
%                       'FLUX_PEAK' - The intensity value in the peak.
%                       'FLUX' - The oveall intensity in the segment. 
%            'IntType' - integer class - used for the X_PEAK,Y_PEAK 
%                        computations. default is iunit16 
%                        (it is highly recommended to use the default IntType 
%                         in order to prevent incompatibility between integers 
%                         types when the number of segments is lower than 255).
%           
%            --- Additional parameters
%            Any additional key,val, that are recognized by one of the
%            following programs:
%            images2sim.m, image2sim.m
% Output : - A structure array of the computed parameters (as indicated in
%            OutPar).  
% Tested : Matlab R2013b
%     By : Tali Engel                      Jan 2015
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: [Moments]=sim_seg_moment(Sim1,Sim_Segment);
% Reliable: 
%--------------------------------------------------------------------------

%FunName = mfilename;

ImageField = 'Im';
SegImField = 'Im';
CatField   = 'Cat';
CatColField = 'Col';
CatColCellField = 'ColCell';

% defaults:
DefV.OutPar = {'Iseg','X','Y','X2','Y2','XY','ALPHA_J2000','DELTA_J2000','NPIX','X_PEAK','Y_PEAK','ALPHA_PEAK','DELTA_PEAK','FLUX_PEAK','FLUX'};
DefV.IntType = @uint16;
DefV.ImageField = ImageField;
DefV.SegImField = SegImField;


InPar = set_varargin_keyval(DefV,'n','use',varargin{:});

ImageField = InPar.ImageField;
SegImField = InPar.SegImField;

%--- read images ---
Sim = images2sim(Sim,varargin{:});
Nim = numel(Sim);

SimSeg = images2sim(SimSeg,varargin{:});
NimSeg = numel(SimSeg);

if ((Nim~=NimSeg) && (Nim~=1) && (NimSeg~=1))
    error('Illegal number of segmenation images');
end

N = max(Nim,NimSeg); % in case one of them equals 1.

%--- define boolean variables according to InPar.OutPar --- 

% 1. First moment (X,Y).
if (any(strcmp('X',InPar.OutPar)))
    X1_out = true;
else
    X1_out = false;
end

if (any(strcmp('Y',InPar.OutPar)))
    Y1_out = true;
else
    Y1_out = false;
end

% 2. Second moment (X2,Y2,XY).
if (any(strcmp('X2',InPar.OutPar)))
    X2_out = true;
else
    X2_out = false;
end

if (any(strcmp('Y2',InPar.OutPar)))
    Y2_out = true;
else
    Y2_out = false;
end

if (any(strcmp('XY',InPar.OutPar)))
    XY_out = true;
else
    XY_out = false;
end

% 3. first moment in celestial coordinates (RA,DEC)
%   (ALPHA_J2000,DELTA_J2000).
if (any(strcmp('ALPHA_J2000',InPar.OutPar)))
    Long1_out = true;
else
    Long1_out = false;
end

if (any(strcmp('DELTA_J2000',InPar.OutPar)))
    Lat1_out = true;
else
    Lat1_out = false;
end

% 4. Number of pixels (NPIX).
if (any(strcmp('NPIX',InPar.OutPar)))
    Npix_out = true;
else
    Npix_out = false;
end   
     
% 5. Position of the pixel with maximum intensity
%    (X_PEAK,Y_PEAK).                            
if (any(strcmp('X_PEAK',InPar.OutPar)))
    XPeak_out = true;
else
    XPeak_out = false;
end

if (any(strcmp('Y_PEAK',InPar.OutPar)))
    YPeak_out = true;
else
    YPeak_out = false;
end

% 6. Position of the pixel with maximum intensity, in celestial coordinates
%    (ALPHA_PEAK,DELTA_PEAK).
         
if (any(strcmp('ALPHA_PEAK',InPar.OutPar)))
    LongP_out = true;
else
    LongP_out = false;
end

if (any(strcmp('DELTA_PEAK',InPar.OutPar)))
    LatP_out = true;
else
    LatP_out = false;
end

% 7. The intensity value in the peak (FLUX_PEAK).
if (any(strcmp('FLUX_PEAK',InPar.OutPar)))
    ValPeak_out = true;
else
    ValPeak_out = false;
end
   
% 8. The oveall intensity of the segment (FLUX).
if (any(strcmp('FLUX',InPar.OutPar))),
    flux_out = true;
else
    flux_out = false;
end


% a loop over all images - compute output parameters for each segment
Npar = numel(InPar.OutPar);

Moments = struct_def(InPar.OutPar,N,1);
for I=1:1:N,
    Iim   = min(I,Nim);           % images counter
    Iseg  = min(I,NimSeg);        % segmentation images counter
        
    % general defenitions
    Im = Sim(Iim).(ImageField);
    Seg_Im = SimSeg(Iseg).(SegImField);
    
    [Ly,Lx] = size(Im);
    [Lys,Lxs] = size(Seg_Im);
    if ((Ly~=Lys) || (Lx~=Lxs))
        error('Sizes of image and segmentation image are not equal, for image number = %d',I);
    end    
    
    [Matx,Maty] = meshgrid(1:Lx,1:Ly);
    % segments #0 - surrounds every segment (-->-1)
    %          #1 - represents the bg       (-->0)
    %          #2,3,... - represents image segments   (-->1,2,...)
    % in order to prevent confusen 1 is substracted from all the segmented
    % image elements:
    Seg_Im = Seg_Im - 1;   % background is now indicated as 0.

    % number of segments
    Numseg = max(Seg_Im(:));      

    Cat(Iim).(CatField)        = zeros(Numseg,Npar);
    Cat(Iim).(CatColCellField) = InPar.OutPar;
    Cat(Iim).(CatColField)     = cell2struct(num2cell(1:1:Npar),InPar.OutPar,2);
    if (any(strcmp('Iseg',InPar.OutPar))),
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).Iseg) = (1:1:Numseg).';
    end
    
    % max number of pixels in a segment
    Max_count = max(histc(Seg_Im(:),1:1:Numseg));

    % sort the pixels by segments (column wise).
    % initialization
    Seg_mat = zeros(Max_count,Numseg);
    Is      = zeros(Max_count,Numseg);
    Xs      = zeros(Max_count,Numseg);
    Ys      = zeros(Max_count,Numseg);
    
    % loop over all pixels
    Counter = zeros(Numseg,1);
    Pind    = find(Seg_Im>0);     %returns a vector containing the linear indices of the 
                                  %segments (the positive elements in the array).
    % build a matrix (Seg_mat) of the linear indeces of the segments'
    % pixels, each segment "members" in a different column.
    for Ipix=1:length(Pind)
        Lpix = Pind(Ipix);               % linear index
        SegVal = Seg_Im(Lpix);           % the segment value 
        Counter(SegVal) = Counter(SegVal)+1; 
        Row = Counter(SegVal);     
        Seg_mat(Row,SegVal) = Lpix;      
    end                                

    % build three matrices: Is, Xs, Ys
    Npix = zeros(1,Numseg);
    for Inumseg=1:1:Numseg
        Col         = Seg_mat(:,Inumseg);
        Len         = sum(Col>0);
        Is(1:Len,Inumseg) = Im(Col(1:Len));
        Xs(1:Len,Inumseg) = Matx(Col(1:Len));
        Ys(1:Len,Inumseg) = Maty(Col(1:Len));
        Npix(Inumseg)     = Len;                % number of pixels in each segment
    end
    
    % 1. First moment (X,Y).
    if (X1_out)
        X1 = sum(Is.*Xs)./sum(Is);
        %Moments(I).X = X1;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).X) = X1;
    end
    if (Y1_out)
        Y1 = sum(Is.*Ys)./sum(Is);
        %Moments(I).Y = Y1;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).Y) = Y1;
    end
    
    % 2. Second moment (X2,Y2,XY).
    if (X2_out)
        %X2 = (sum(Is.*Xs.^2)./sum(Is))-X1.^2; 
        %Moments(I).X2 = X2;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).X2) = (sum(Is.*Xs.^2)./sum(Is))-X1.^2; 
    end
    if (Y2_out)
        %Y2 = (sum(Is.*Ys.^2)./sum(Is))-Y1.^2;
        %Moments(I).Y2 = Y2;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).Y2) = (sum(Is.*Ys.^2)./sum(Is))-Y1.^2;
    end
    if (XY_out)
        %XY = (sum(Is.*Xs.*Ys)./sum(Is))-X1.*Y1;
        %Moments(I).XY = XY;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).XY) = (sum(Is.*Xs.*Ys)./sum(Is))-X1.*Y1;
    end
    
    % 3. first moment in celestial coordinates (RA,DEC)
    %    (ALPHA_J2000,DELTA_J2000).
    if ( Long1_out || Lat1_out );
        % convert X1,Y1 to [RA DEC]
        [Long1,Lat1] = xy2sky(Sim(Iim),X1,Y1);
        if (Long1_out)
            %Moments(I).ALPHA_J2000 = Long1';
            Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).ALPHA_J2000) = Long1.';
        end
        if (Lat1_out)
            %Moments(I).DELTA_J2000 = Lat1';
            Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).DELTA_J2000) = Lat1.';
        end
    end    
    
    % 4. Number of pixels (NPIX).
    if (Npix_out)
        Moments(I).NPIX = Npix;
        Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).NPIX) = Npix;
    end
    
    % Peak with maximum intensity
    if (XPeak_out || YPeak_out || LongP_out || LatP_out || ValPeak_out)
        [ValPeak,IndPeak] = max(Is);   
        
        % 5. Position of the pixel with maximum intensity
        %    (X_PEAK,Y_PEAK).                            
        if (XPeak_out)
%            XPeak = Xs(sub2ind(uint16(size(Is)),uint16(IndPeak),uint16(1:Numseg)));
            XPeak = Xs(sub2ind(InPar.IntType(size(Is)),InPar.IntType(IndPeak),InPar.IntType(1:Numseg)));            
            %Moments(I).X_PEAK = XPeak;
            Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).X_PEAK) = XPeak;
        end
        if (YPeak_out)
%            YPeak = Ys(sub2ind(uint16(size(Is)),uint16(IndPeak),uint16(1:Numseg)))
            YPeak = Ys(sub2ind(InPar.IntType(size(Is)),InPar.IntType(IndPeak),InPar.IntType(1:Numseg)));
            %Moments(I).Y_PEAK = YPeak;
            Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).Y_PEAK) = YPeak;
        end
         
        % 6. Position of the pixel with maximum intensity, in celestial coordinates
        %    (ALPHA_PEAK,DELTA_PEAK).
        if (LongP_out || LatP_out),
            [LongP,LatP]=xy2sky(Sim(Iim),XPeak,YPeak);
            if (LongP_out)
                %Moments(I).ALPHA_PEAK = LongP';
                Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).ALPHA_PEAK) = LongP.';
            end
            if (LatP_out)
                %Moments(I).DELTA_PEAK = LatP';
                Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).DELTA_PEAK) = LatP.';
            end
        end
         
        % 7. The intensity value in the peak (FLUX_PEAK).
        if (ValPeak_out)
            Moments(I).FLUX_PEAK = ValPeak;
            Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).FLUX_PEAK) = ValPeak;
        end
    end
     
     % 8. The oveall intensity of the segment (FLUX).
     if (flux_out)
         %Moments(I).FLUX = sum(Is);
         Cat(Iim).(CatField)(:,Cat(Iim).(CatColField).FLUX) = sum(Is);
     end
     
end

    
     
     
     
     
     
     
