Speeding up Image Processing using Parallel Computing

Hi all - I'm segmenting large (16k-20k x 16k-20k pixels) histology images based on the Euclidean distance to numerous known colors. In my current application, I'm segmenting various components of healing bone from slides stained with Goldner's Trichrome. I'm currently doing this in a double for loop. As you can imagine, this is very process intensive, and some of my samples are running 12-15 minutes. I have the parallel computing toolbox and a good multi-core computer, so I'd like to start implementing parallel computing into my image analysis. When I make any of these loops a "parfor" loop, I get errors that parfor cannot work in this type of loop. Could anybody help me out as to how I can speed up my application with parallel computing? Thanks, Tristan

A=imread('sample.tif');
cart=double(A);
% Assign colors of interest in arrays 
minBoneRGB=[50.54, 169.34, 169.17]; 
unminBoneNormalRGB=[168, 71.5, 122.2]; 
unminBoneBrightRGB=[173.3, 23.3, 89.8]; 
immatureBoneRGB = [130.37, 179, 172.34]; 
backgroundRGB=[204.37, 199.5, 190.6];

    minBone=zeros(size(cart,1), size(cart,2));
    unminBoneNormal=zeros(size(cart,1), size(cart,2));
    unminBoneBright=zeros(size(cart,1), size(cart,2));
    immatureBone=zeros(size(cart,1), size(cart,2));
    background=zeros(size(cart,1), size(cart,2));

    % Determine intensity of each tissue in each channel
    for i=1:size(cart,1)
        for j=1:size(cart,2)
            if VOL2(i,j)==1
                minBone(i,j)=255-sqrt(((minBoneRGB(1)-cart(i,j,1))^2)+((minBoneRGB(2)-cart(i,j,2))^2)+((minBoneRGB(3)-cart(i,j,3))^2));
                if minBone(i,j)<0
                    minBone(i,j)=0;
                end
                  unminBoneNormal(i,j)=255-sqrt(((unminBoneNormalRGB(1)-cart(i,j,1))^2)+((unminBoneNormalRGB(2)-cart(i,j,2))^2)+((unminBoneNormalRGB(3)-cart(i,j,3))^2));
                  if unminBoneNormal(i,j)<0
                      unminBoneNormal(i,j)=0;
                  end
                   unminBoneBright(i,j)=255-sqrt(((unminBoneBrightRGB(1)-cart(i,j,1))^2)+((unminBoneBrightRGB(2)-cart(i,j,2))^2)+((unminBoneBrightRGB(3)-cart(i,j,3))^2));
                  if unminBoneBright(i,j)<0
                      unminBoneBright(i,j)=0;
                  end
                  immatureBone(i,j)=255-sqrt(((immatureBoneRGB(1)-cart(i,j,1))^2)+((immatureBoneRGB(2)-cart(i,j,2))^2)+((immatureBoneRGB(3)-cart(i,j,3))^2));
                  if immatureBone(i,j)<0
                      immatureBone(i,j)=0;
                  end
                  background(i,j)=255-sqrt(((backgroundRGB(1)-cart(i,j,1))^2)+((backgroundRGB(2)-cart(i,j,2))^2)+((backgroundRGB(3)-cart(i,j,3))^2));
                  if background(i,j)<0
                      background(i,j)=0;
                  end
              end
          end
      end
      end

Answer 1

[Edric Ellis]

I think vectorization using BSXFUN is likely to help you out more here than parallelization if I've understood your problem correctly. Here's my stab at doing that, in a simplified sort of way, but hopefully you can apply it to your problem (shouldn't be too far away).

    % Set up a sample 'image'
    img = rand(5, 5, 3); img(3, 3, :) = 0.5

    % Set up the thing we compare against
    compare = [0.5, 0.5, 0.5];

    % Ensure 'compare' has the right shape for the impending BSXFUN
    compare = reshape(compare, 1, 1, 3);

    % Subtract compare from img, expanding compare as needed
    simpleDiff = bsxfun(@minus, img, compare);

    % Compute the resulting distance
    distance = sum( simpleDiff .* simpleDiff, 3 )

I've missed the clamping to zero piece, but that should be simple doing something like

    distance(distance<0) = 0;

Answer 2

[Tristan]

Edric -

thanks for your input. This is an interesting and very fast way of getting the same relative relationship of the "distance" between the color of interest and each pixel in an image. I certainly think I can implement it in my analyses. If you have a few more minutes, could you please explain the fifth line of your syntax where distance is computed? I'm not fully grasping how the sum of the squares of the third dimension of simpleDiff gets us to a variation of the Euclidean distance between two RGB colors. I see that it's working wonderfully but I'd like to fully understand it. Thanks,

Tristan

[Edric Ellis]

Hi Tristan, the BSXFUN call creating simpleDiff expands the array 'compare' so that

simpleDiff(i,j,k) = img(i,j,k) - compare(k)

for all i,j,k. Summing the elementwise square of simpleDiff along the third dimension means that

distance(i,j) = simpleDiff(i,j,1) .* simpleDiff(i,j,1) + ...
  simpleDiff(i,j,2) .* simpleDiff(i,j,2) + ...
  simpleDiff(i,j,3) .* simpleDiff(i,j,3);

Answer 3

[Jan]

I'd start with simplifying the code at first. E.g. size(cart,1) is a function, which consumes a tiny, but measureable piece of time. So better define this once outside the loops and store it in a temporary variable. As side-effect ths code will become much leaner and locations of further simplifications will appear.