Asked by Tristan
on 21 Sep 2012

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 by Edric Ellis
on 21 Sep 2012

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 by Tristan
on 21 Sep 2012

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
on 24 Sep 2012

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 by Jan Simon
on 24 Sep 2012

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.

Opportunities for recent engineering grads.

## 0 Comments