Обратите внимание, что причина того, что вам не хватает памяти, вероятно, в том, что вы использовали parfor для замены одного из двух циклов в вашем коде. В обоих случаях каждый рабочий поток будет создавать массив того же размера, что и foto во время обработки. Обратите внимание, что во внутреннем цикле sum(L(:) == f) создает логический массив размером L, а затем суммирует его значения (я не думаю, что JIT достаточно умен, чтобы оптимизировать этот промежуточный массив).
Короче говоря, распараллеливание операции над таким большим изображением, как вы сделали, нецелесообразно. Правильный способ его распараллеливания состоит в том, чтобы разрезать изображение на плитки и обрабатывать каждую плитку в отдельном потоке. Если фрагменты маленькие (предположение, которое я осмелюсь сделать, учитывая имя), то должна быть возможность обрабатывать плитки, используя только небольшое перекрытие (плитки должны перекрываться так, чтобы каждый фрагмент был полностью внутри, по крайней мере, на плитке). В этом сценарии немного сложнее удалить дубликаты, поэтому распараллеливание не является тривиальным. Однако я надеюсь, что приведенное ниже предложение вообще не будет необходимости распараллеливать код.
Из кода видно, что фрагменты одного и того же класса не соприкасаются. Но не ясно, что фрагменты из разных классов не соприкасаются (т. Е. Код выдает одинаковый результат, если они это сделали). При условии, что они не , можно избежать обеих петель.
Идея состоит в том, чтобы пометить все фрагменты, независимо от класса, за один раз. Поэтому вместо того, чтобы вызывать bwlabeln один раз для каждого класса, вы вызываете его только один раз. Я не знаю, сколько существует классов, но это потенциально значительное сокращение времени вычислений.
Далее используйте regionprops, чтобы определить для каждого фрагмента его размер и класс. Эта операция в принципе также может быть выполнена путем итерации по изображению только один раз. Обратите внимание, что ваш код FragmentSize(f,1) = sum(L(:) == f) перебирает изображение один раз для каждого фрагмента. Учитывая размер изображения, могут быть миллионы фрагментов. Это может обеспечить сокращение времени на 6 порядков.
С этого момента мы имеем дело только с выводом regionprops, который может содержать (на порядок) миллион элементов, тривиальную величину (на 3 порядка меньше, чем количество пикселей).
Это может быть код:
L = bwlabeln(foto>0);
cls = regionprops(L,foto,'MinIntensity','Area');
clear L
sz = [cls.Area];
cls = [cls.MinIntensity];
NumberOfClasses = max(cls);
FRAGMENTSIZESCLASS = cell(NumberOfClasses,1);
for ii=1:NumberOfClasses
FRAGMENTSIZESCLASS{ii} = sz(cls==ii);
end
Этот последний цикл может не понадобиться, я не нашел быстрой альтернативы. Я не могу себе представить, что это дорого, но если это так, тривиально распараллелить или улучшить его, отсортировав cls и используя diff, чтобы найти индексы, с которых начинается новый класс.
Можно переписать приведенный выше код, предложив @ bla bwconncomp. Эта функция возвращает структуру, содержащую массив ячеек с индексами для всех пикселей с каждой меткой. Тогда нет необходимости использовать regionprops, можно непосредственно найти размер (как показало @bla) и использовать первый индекс для каждой метки, чтобы найти класс (путем индексации в foto):
cc = bwconncomp(foto>0);
sz = cellfun(@numel,cc.PixelIdxList);
cls = cellfun(@(indx)foto(indx(1)),cc.PixelIdxList);
NumberOfClasses = max(cls);
FRAGMENTSIZESCLASS2 = cell(NumberOfClasses,1);
for ii=1:NumberOfClasses
FRAGMENTSIZESCLASS2{ii} = sz(cls==ii);
end
Это было в 3-4 раза быстрее для небольшого тестового изображения 256x256 с 63 фрагментами. Однако, учитывая размер изображения, с которым вы имеете дело, я боюсь, что это может быть очень неэффективно. Единственный способ узнать это - попробовать оба подхода и рассчитать время!
Несколько замечаний о вашем коде:
FRAGMENTSIZESCLASS = struct([]);
Вы инициализируете его как пустой структурный массив, но затем используете {} для индексации в нем, преобразовывая его в массив ячеек. Всегда хорошо предварительно распределить массивы, как я делал выше:
FRAGMENTSIZESCLASS = cell(NumberOfClasses,1);
NumberFragments=max(max(L));
Это создает максимальную проекцию L на горизонтальную ось (80k элементов), а затем находит максимум в этом. Более эффективно изменить форму матрицы, как вы это делали в другом месте:
NumberFragments = max(L(:));