Применение симплекс-метода

function [Xr, I1, MeanImg] = simplex(Xin, J, I)
% Xin – анализируемые факторы;
% J   – интервал варьирования факторов;
% I   – изображение.
%% Подготовка матрицы
k = size( Xin, 1 );
p = ( k - 1 + sqrt( k + 1 )) / ( k * sqrt( 2 ));
q = ( sqrt( k + 1 ) - 1 ) / ( k * sqrt( 2 ));
X = zeros( k ) + q;
for i = 1 : k, X( i, i ) = p; end
X = [ zeros(1, k ); X ];
X = repmat( transp( Xin ), k + 1, 1 ) + X .* repmat( transp( J ), k + 1, 1);
%% Эксперимент
S = zeros( k + 1, 1 );
for i = 1 : k + 1
    T = double( uint8( X( i, : )' ));
    I1 = uint8( rangfilt( I, T ));
    S( i ) = contrastBC( I1 );
end
Xr = [ X S ];
rp = 4; % счетчик (номер текущей новой вершины)
while rp < 50
    sn = find( S == min( S ), 1 ); % sn – номер строки с наихудшим откликом
    for p = 1 : k + 1
        if sum( sort( X( p, :))~=X( p, : )) > 0
            sn = p;
            break;
        end
    end

    X( sn, : ) = ( 2 / k ) * sum( X ) - ( 2 / k + 1 ) * X( sn, : );
    T = double( uint8( X( sn, : )' ));
    I1 = uint8( rangfilt( I, T ));
    [S( sn ), MeanImg] = contrastBC( I1 );
    Xr = [ Xr; [ X( sn, : ) S( sn ) ]]; % Xr – Матрица, хранящая значения каждого шага
    rp = rp + 1;
end