function Y = apply_mask(X, X_axes, varargin)
% Y = apply_mask(X, X_axes, option/value pairs)
%   apply mask (as indicated with option/value pairs) to spectrum X (with X
%   and X_axes as read from an NMR Pipe format file using read_nmrp).
%   Option / value pairs include:
%   
%   'water' [w1 w2]        mask water line in direct (first) dimension between chemical shifts w1 and w2
%   'donor' mask_array     interpolate mask_array to apply to each donor slice and apply to each donor plane  
%   'acceptor' mask_array  interpolate mask_array to apply to each acceptor slice and apply to each acceptor plane
%   'operation' op         op is a function handle such that op(A,B) = C for arrays of arbitrary dim. and size
%                          the masked spectrum is defined as Y = op(X,mask) and the default is min
%                          note: ideally operation(operation(a,b),c) = operation(operation(a,c),b)  
%   'regularize' alpha     add alpha to the diagonal elements of the matrix to ensure the 'target' spectrum to be 
%                          reconstructed via covariance is positive-semidefinite
%   'footprint'  s2n       foot-printing (specialized acceptor mask, therefore, this option cannot be used with the 
%                         'acceptor' option) with critical value of s2n above the median of the absolute value of the
%                         noise floor

% set defaults

operation = @(a,b)(sign(a).*min(abs(a),b)); % assume mask is non-negative
water = 0;
donor = 0;
acceptor = 0;
alpha = 0;

num_args_in = nargin;
if(num_args_in < 2)
    help apply_mask
    Y = 0;
    return
end

cur_dim_info = squeeze(X_axes(27,:)); % vector with 1 for acceptor dims
cur_size_row = squeeze(X_axes(26,:));
cur_acceptor_dims = find(cur_dim_info == 1);
cur_donor_dims = find(cur_dim_info == -1);

num_dims = sum(cur_size_row > 1);


for(i=1:2:(num_args_in-2))
    opt_name = varargin{i};
	opt_value = varargin{i+1};
    
	switch opt_name(1) % first char sufficient to ID option type
		case {'w','W'} % water line to remove
            water = opt_value;
		case {'d','D'} % mask array to apply to donor slices 
			donor = opt_value;
        case {'A','a'} % mask array to apply to donor slices 
			acceptor = opt_value; 
		case {'o','O'} % handle to function containing operation to use
			operation = opt_value;
        case {'r','R'}
            alpha = opt_value;
        case {'f','F'}
            crit_s2n = opt_value;
            acceptor = abs(X);

            for(i=length(cur_donor_dims):-1:1)
                cur_don_dim = cur_donor_dims(i);
                acceptor = squeeze(max(acceptor,[],cur_don_dim));
            end

            noise_level = median(acceptor(:));
            acceptor = acceptor > (noise_level.*crit_s2n);
             
        % add more cases here as options are expanded
    end
end

Y = X;
clear X;

if(alpha ~= 0) % regularize spectrum
    axes_ranges = ...
        [(X_axes(16,:)./X_axes(14,:))' (X_axes(16,:)./X_axes(14,:)+X_axes(20,:)./X_axes(14,:))' (X_axes(20,:)./(X_axes(14,:).*cur_size_row))'];
       %           min                                 max                                          inc / spacing
    pair_idx = 0;

	for(i=1:length(cur_donor_dims))
		for(j=1:length(cur_acceptor_dims))
			donor_idx = cur_donor_dims(i);
			accpt_idx = cur_acceptor_dims(j);
			if((axes_ranges(donor_idx,1) < axes_ranges(accpt_idx,2))&&(axes_ranges(accpt_idx,1) < axes_ranges(donor_idx,2)))
				pair_idx = pair_idx + 1;
				da_pairs(pair_idx,:) = [donor_idx accpt_idx];			
			end
		end
	end   

	for(i=1:4)
		ranges2reg{i} = 1:cur_size_row(i);
    end
    
    a1_icds = mod(round(...
			(((1:cur_size_row(da_pairs(1,1))).*axes_ranges(da_pairs(1,1),3)) + axes_ranges(da_pairs(1,2),1) - axes_ranges(da_pairs(1,1),1))/...
				axes_ranges(da_pairs(1,2),3))-1,cur_size_row(da_pairs(1,2)))+1;
    if(size(da_pairs,1)>1) % max size of pairs list should be 2 then
        a2_icds = mod(round(...
			(((1:cur_size_row(da_pairs(2,1))).*axes_ranges(da_pairs(2,1),3)) + axes_ranges(da_pairs(2,2),1) - axes_ranges(da_pairs(2,1),1))/...
				axes_ranges(da_pairs(2,2),3))-1,cur_size_row(da_pairs(2,2)))+1;
    end

    for(i=1:cur_size_row(da_pairs(1,1)))
		ranges2reg{da_pairs(1,1)} = i;
        ranges2reg{da_pairs(1,2)} = a1_icds(i);
        
        if(size(da_pairs,1) > 1) % max size of pairs list should be 2, so only two cases to consider
			for(j=1:cur_size_row(da_pairs(2,1)))

                ranges2reg{da_pairs(2,1)} = j;
				ranges2reg{da_pairs(2,2)} = a2_idcs(j);

				Y(ranges2reg{1},ranges2reg{2},ranges2reg{3},ranges2reg{4}) = ...
					Y(ranges2reg{1},ranges2reg{2},ranges2reg{3},ranges2reg{4}) + alpha;					
			end
        else
			Y(ranges2reg{1},ranges2reg{2},ranges2reg{3},ranges2reg{4}) = ...
				Y(ranges2reg{1},ranges2reg{2},ranges2reg{3},ranges2reg{4}) + alpha;
		end
    end
end 

if(numel(water) > 1)
    water = ppm2inc([water(1);water(2)], X_axes);
    water = sort(water);
    Y(water(1):water(2),:) = operation(Y(water(1):water(2),:), 0); % don't need to worry about dimensionality of
end                                                                % Y, this routine works regardless 

if(numel(donor) > 1) %% note that 3D donor/acceptor masks not fully tested (as of March 2010) in the context of an actual spec masking
    
    Y_permuted = permute(Y, [cur_acceptor_dims cur_donor_dims]); % place acceptor dims first
    clear Y; % save space
	Y_reshaped = reshape(Y_permuted,[prod(cur_size_row(cur_acceptor_dims)) prod(cur_size_row(cur_donor_dims))]);  
                                                     % reshape array into
                                                     % an acceptors x donors matrix
    clear Y_permuted; % save space
    
	switch(sum(size(donor) > 1)) % number of donor dimensions
		case {1}
			interpX1 = linspace(1,max(size(donor)),cur_size_row(cur_donor_dims(1)));
			donor = interp1(donor,interpX1);
			donor = donor(:)'; 
		case {2}
			interpX1 = linspace(1,size(donor,1),cur_size_row(cur_donor_dims(1)));
			interpX2 = linspace(1,size(donor,2),cur_size_row(cur_donor_dims(2)));
			donor = interp2(donor,interpX1,interpX2'); 
			donor = donor(:)'; 
    	case {3} 
			interpX1 = linspace(1,size(donor,1),cur_size_row(cur_donor_dims(1)));
			interpX2 = linspace(1,size(donor,2),cur_size_row(cur_donor_dims(2)));
			interpX3 = linspace(1,size(donor,3),cur_size_row(cur_donor_dims(3)));
            [X1 X2 X3] = meshgrid(interpX1,interpX2,interpX3);
			donor = interp3(donor,X1,X2,X3);
			donor = donor(:)'; 
    end
    
    N_acc = size(Y_reshaped,1);
    for(i=1:N_acc)
        Y_reshaped(i,:) = operation(Y_reshaped(i,:),donor);
    end
    
    Y_permuted = reshape(Y_reshaped,[cur_size_row(cur_acceptor_dims) cur_size_row(cur_donor_dims)]);
    clear Y_reshaped;
    Y = ipermute(Y_permuted, [cur_acceptor_dims cur_donor_dims]); % reverse everything
    clear Y_permuted;
end

if(numel(acceptor) > 1)

    Y_permuted = permute(Y, [cur_acceptor_dims cur_donor_dims]); % place acceptor dims first 
    clear Y; % save space
	Y_reshaped = reshape(Y_permuted,[prod(cur_size_row(cur_acceptor_dims)) prod(cur_size_row(cur_donor_dims))]);  
                                                     % reshape array into
                                                     % an acceptors x donors matrix
    clear Y_permuted; % save space
    
    switch(sum(size(acceptor) > 1)) % number of acceptor dimensions
		case {1}
			interpX1 = linspace(1,max(size(acceptor)),cur_size_row(cur_acceptor_dims(1)));
			acceptor = interp1(acceptor,interpX1);
			acceptor = acceptor(:); 
		case {2}
			interpX1 = linspace(1,size(acceptor,1),cur_size_row(cur_acceptor_dims(1)));
			interpX2 = linspace(1,size(acceptor,2),cur_size_row(cur_acceptor_dims(2)));
			acceptor = interp2(acceptor,interpX1,interpX2'); 
			acceptor = acceptor(:); 
    	case {3} 
			interpX1 = linspace(1,size(acceptor,1),cur_size_row(cur_acceptor_dims(1)));
			interpX2 = linspace(1,size(acceptor,2),cur_size_row(cur_acceptor_dims(2)));
			interpX3 = linspace(1,size(acceptor,3),cur_size_row(cur_acceptor_dims(3)));
            [X1 X2 X3] = meshgrid(interpX1,interpX2,interpX3);
			acceptor = interp3(acceptor,X1,X2,X3);
			acceptor = acceptor(:); 
    end
    
    N_don = size(Y_reshaped,2);
    for(j=1:N_don)
        Y_reshaped(:,j) = operation(Y_reshaped(:,j),acceptor);
    end
    
    Y_permuted = reshape(Y_reshaped,[cur_size_row(cur_acceptor_dims) cur_size_row(cur_donor_dims)]);
    clear Y_reshaped;
    Y = ipermute(Y_permuted, [cur_acceptor_dims cur_donor_dims]); % reverse everything
    clear Y_permuted;
end