MATLAB: Need to accelerate an operation on a matrix

digital image processingMATLABparallel computingParallel Computing Toolboxspeed

I'm writing code that involves looping through each element in arrays with appx. 10e5 elements. For each value in the array, a filter is applied to the cells around it. This filtered neighborhood is used to determine the value assigned to a value in a new matrix. This process is the bottleneck in the code, so I'd appreciate any help in making this loop faster. I tried a gpuArray-based implementation but it was significantly slower than what I'm currently doing. One thing I find curious is that a simple indexing of the large array takes longer than some of the more complex operations.
Here's the code:
for k = 1:ncells
    pos = allinds(k,:);
    yf = pos(1) + sensdist;
    xf = pos(2) + sensdist;
    % Slow operation: apply the neighborhood filter
    neighbors = subfield(yf-sensdist:yf+sensdist, xf-sensdist:xf+sensdist) .*obj.nhood;
    %Slow operation: look up next value for this index
    nextstate = trans(1,subfield(yf,xf));
    %Slow operation: check how many neighbors possess are the next value
    quorum = sum(sum(neighbors == nextstate));
    if quorum > 1;
      if ~isempty(zipintersect(quorum,obj.go))
        temp(k) = nextstate;
      end
    end
  end

Best Answer

This may help some. The key points are,
  • Pre-allocate temp if you're not doing so already.
  • Create a table of NextStates in advance for all positions (yf,xf). Makes it simpler to access
  • Get rid of 2D indexing and replaces it by linear indexing.
  • Get rid of double sum() by shaping the neighborhood as a vector all the time instead of a 2D submatrix.
  • Get rid of double "if"
    sz=size(subfield);
    w=sensdist+1;
    ww=w+sensdist;
    [ii,jj]=ndgrid(1:ww);
    jumps= sub2ind(sz,ii,jj) - sub2ind(sz, w, w) ;
      jumps=jumps(:);
      obj.nhood=obj.nhood(:);  
    NextStates=reshape( trans(1,subfield(:)), sz ) ;
    allpos=sub2ind(sz,allinds(:,1)+sensdist,allinds(:,2)+sensdist);
    temp=nan(1,ncells);
    for k = 1:ncells
        pos=allpos(k);
        % Slow operation: apply the neighborhood filter
        neighbors = subfield(pos+jumps) .*obj.nhood;
        %Slow operation: look up next value for this index
        nextstate = NextStates(pos);
        %Slow operation: check how many neighbors possess are the next value
        quorum = sum(neighbors == nextstate);
        if quorum > 1 && ~isempty(zipintersect(quorum,obj.go))
            temp(k) = nextstate;
        end
      end
Related Question