Обработка изображений: алгоритм занимает слишком много времени в MATLAB - PullRequest
Новая
       53

Обработка изображений: алгоритм занимает слишком много времени в MATLAB

0 голосов
/ 29 мая 2018

Я работаю в MATLAB для обработки двух изображений 512x512, изображения домена и изображения диапазона.Я пытаюсь сделать следующее:

  1. Разделить изображения домена и диапазона на блоки 8x8 пикселей
  2. Для каждого блока 8x8 в изображении домена я должен применить линейныйпреобразования в него и сравнить каждый из 4096 преобразованных блоков с каждым из 4096 блоков диапазона.
  3. Вычислить ошибку в каждом случае между преобразованным блоком и блоком изображения диапазона и найти минимальную ошибку.
  4. Наконец, у меня будет для каждого блока диапазона 8x8 идентификатор блока домена 8x8для которого ошибка была минимальной (ошибка между блоком диапазона и преобразованным блоком домена)

Для достижения этого я написал следующий код: 

RangeImagecolor = imread('input.png'); %input is 512x512
DomainImagecolor = imread('input.png'); %Range and Domain images are identical

RangeImagetemp = rgb2gray(RangeImagecolor);
DomainImagetemp = rgb2gray(DomainImagecolor);

RangeImage = im2double(RangeImagetemp);

DomainImage = im2double(DomainImagetemp);
%For the (k,l)th 8x8 range image block
for k = 1:64
    for l = 1:64
        minerror = 9999;
        min_i = 0;
        min_j = 0;

        for i = 1:64
            for j = 1:64

                 %here I compute for the (i,j)th domain block, the transformed domain block stored in D_trans

                 error = 0;
                 D_trans = zeros(8,8);
                 R = zeros(8,8); %Contains the pixel values of the (k,l)th range block
                 for m = 1:8
                     for n = 1:8
                         R(m,n) = RangeImage(8*k-8+m,8*l-8+n); 
                         %ApplyTransformation can depend on (k,l) so I can't compute the transformation outside the k,l loop.
                         [m_dash,n_dash] = ApplyTransformation(8*i-8+m,8*j-8+n);
                         D_trans(m,n) = DomainImage(m_dash,n_dash);
                         error = error + (R(m,n)-D_trans(m,n))^2; 
                     end 
                 end
                 if(error < minerror)
                     minerror = error; 
                     min_i = i;
                     min_j = j;
                 end
            end
        end 
    end
end

В качестве примераПрименив преобразование, можно использовать преобразование идентичности: 

 function [x_dash,y_dash] = Iden(x,y)
    x_dash = x;
    y_dash = y;
end

Теперь проблема, с которой я сталкиваюсь, - это большое время вычислений.Порядок вычислений в приведенном выше коде составляет 64 ^ 5, что порядка 10 ^ 9.Это вычисление должно занять в худшие минуты или час.На вычисление всего 50 итераций уходит около 40 минут.Я не знаю, почему код работает так медленно.

Спасибо, что прочитали мой вопрос.

Ответы [ 2 ]

0 голосов
/ 29 мая 2018

Вы можете использовать im2col* для преобразования изображения в формат столбца, чтобы каждый блок образовывал столбец матрицы [64 * 4096].Затем примените преобразование к каждому столбцу и используйте bsxfun для векторизации вычисления ошибки. 

DomainImage=rand(512);
RangeImage=rand(512);
DomainImage_col = im2col(DomainImage,[8 8],'distinct');
R =  im2col(RangeImage,[8 8],'distinct');
[x y]=ndgrid(1:8);

function [x_dash, y_dash] = ApplyTransformation(x,y)
    x_dash = x;
    y_dash = y;
end

[x_dash, y_dash] = ApplyTransformation(x,y);
idx = sub2ind([8 8],x_dash, y_dash);

D_trans = DomainImage_col(idx,:);       %transformation is reduced to matrix indexing
Error = 0;
for mn = 1:64
    Error = Error + bsxfun(@minus,R(mn,:),D_trans(mn,:).').^2;
end
[minerror ,min_ij]= min(Error,[],2);         % linear index of minimum of each block;
[min_i  min_j]=ind2sub([64 64],min_ij); % convert linear index to subscript

Объяснение :

Наша цель состоит в том, чтобымаксимально сократить количество петель.Для этого нам следует избегать матричной индексации и вместо этого мы должны использовать векторизацию.Вложенные циклы должны быть преобразованы в один цикл.В качестве первого шага мы можем создать более оптимизированный цикл, как здесь: 

min_ij = zeros(4096,1);

for kl = 1:4096                 %%% => 1:size(D_trans,2)
    minerror = 9999;
    min_ij(kl) = 0;
    for ij = 1:4096             %%% => 1:size(R,2)
        Error = 0;
        for mn = 1:64
            Error = Error + (R(mn,kl) - D_trans(mn,ij)).^2;
        end
        if(Error < minerror)
            minerror = Error; 
            min_ij(kl) = ij;
        end
    end
end

Мы можем переставить циклы, и мы можем сделать самый внутренний цикл в качестве внешнего цикла и отделить вычисление минимума отвычисление ошибки. 

% Computation of the error
Error = zeros(4096,4096);
for mn = 1:64
    for kl = 1:4096
        for ij = 1:4096
            Error(kl,ij) = Error(kl,ij) + (R(mn,kl) - D_trans(mn,ij)).^2;
        end
    end
end

% Computation of the min
min_ij = zeros(4096,1);
for kl = 1:4096
    minerror = 9999;
    min_ij(kl) = 0;
    for ij = 1:4096
        if(Error(kl,ij) < minerror)
            minerror = Error(kl,ij); 
            min_ij(kl) = ij;
        end
    end
end

Теперь код упорядочен так, что его лучше всего векторизовать: 

Error = 0;
for mn = 1:64
    Error = Error + bsxfun(@minus,R(mn,:),D_trans(mn,:).').^2;
end

[minerror ,min_ij] = min(Error, [], 2);   
[min_i  ,min_j] = ind2sub([64 64], min_ij);

* Если у вас нет набора инструментов для обработки изображений, болееэффективную реализацию im2col можно найти здесь .

* Весь расчет занимает менее минуты.

0 голосов
/ 29 мая 2018

Перво-наперво - ваш код ничего не делает.Но вы, вероятно, что-то делаете с этим минимумом ошибок и забыли вставить это здесь, или вам все еще нужно кодировать этот бит.Не берите в голову на данный момент.

Одна большая проблема с вашим кодом заключается в том, что вы вычисляете преобразование для блоков размером 64x64 результирующего изображения и исходного изображения.64 ^ 5 итераций сложной операции должны быть медленными.Скорее, вы должны рассчитать все преобразования сразу и сохранить их. 

allTransMats = cell(64);
for i = 1 : 64 
   for j = 1 : 64
      allTransMats{i,j} = getTransformation(DomainImage, i, j)
   end
end
function D_trans = getTransformation(DomainImage, i,j)
D_trans = zeros(8);
for m = 1 : 8
   for n = 1 : 8
      [m_dash,n_dash] = ApplyTransformation(8*i-8+m,8*j-8+n);
      D_trans(m,n) = DomainImage(m_dash,n_dash);
   end
end
end

Это служит для получения allTransMat и находится вне цикла k, l.Желательно в качестве простой функции.

Теперь вы делаете большой цикл k, l, i, j, где сравниваете все элементы по мере необходимости.Сравнение также можно проводить по блокам вместо заполнения небольшой матрицы 8x8, но по какой-то причине делать это для каждого элемента. 

m = 1 : 8;
n = m;
for ...
    R = RangeImage(...); % This will give 8x8 output as n and m are vectors.
    D = allTransMats{i,j};
    difference = sum(sum((R-D).^2));
    if (difference < minDifference) ...
end

Даже если это простой случай без преобразований, это значительно ускоряет код.

Наконец, вы уверены, что вам нужно сравнить каждый блок преобразованного выхода с каждым блоком в источнике?Обычно вы сравниваете block1(a,b) с block2(a,b) - блоками (или пикселями) в одной и той же позиции.

РЕДАКТИРОВАТЬ: для allTransMats требуются также k и l.Уч.НЕТ ПУТИ, чтобы сделать это быстро за одну итерацию, так как вам требуется 64 ^ 5 вызовов ApplyTransformation (или векторизация этой функции, но даже тогда это может быть не быстро - мы должны увидеть эту функцию, чтобы помочь здесь),Поэтому я повторю свой совет, чтобы сгенерировать все преобразования, а затем выполнить поиск: эту верхнюю часть ответа с генерацией allTransMats следует изменить, чтобы иметь все 4 цикла и сгенерировать allTransMats{i,j,k,l};.Это будет медленным, нет способа обойти это, как я упоминал в верхней части правки.Но вы платите один раз, так как после сохранения allTransMats все последующие анализы изображения смогут просто загрузить его, а не генерировать снова.

Но ... что вы даже делаете?Преобразование, которое зависит от индексов блоков источника и назначения, а также индексов пикселей (= всего 6 значений), звучит где-то как ошибка или как основной кандидат на оптимизацию вместо всех остальных.

...