Weird Behaviour of Parallel Processing when Using BlockProc
Show older comments
Hello, I a trying to speed up some analysis of an image (4700 v 128 pixels) that needs breaking up into smaller regions (20x128 pixels) and for each region having a calculation done such as the standard deviation (std2)
I thought the combination of Blockproc and the parallel toolbox would be ideal for this.
However, when I run my analysis, I sometmes get a fast completion time of 0.25s, but other times I run it it can take upto 3s.
1: When I start up my matlab session, i start a parallel pool:
% For parallel Processing (init parpool here)
delete(gcp('nocreate')); % but delete any that are already running
c=parcluster('local');
numwks=c.NumWorkers;
ReportMessage(app,['Parallel processing: NumWorkers = ',num2str(numwks)]);
parpool('Processes',numwks);
p=gcp;
conn=p.Connected;
ReportMessage(app,['...Connected ',num2str(conn)]);
2: I then run my fucntion from a push button callback.
My function first gets an image off a uiaxes, crops it and then performs the blockproc.
% Fast Analysis (SD)
tstart=tic;
ax=app.UIAxes;
IM=getimage(ax);
[sy,sx]=size(IM);
xl=5950; xr=10680; yb=0; yt=0; dx=xr-xl; % This is for image cropping
rect=[xl, yb, xr-xl,sy-yt-yb]; %[x, y, width, height]
IM2=imcrop(IM,rect);
IM=IM2-min(IM2(:)); % Background subtract
[sy,sx]=size(IM);
ColSize=20; % Number of pixels in each ROI - 20 default
nCols=floor(sx/ColSize);
RowSize = 128;
%BlockProc with Parallel
bss = [RowSize,ColSize]; % Each blockproc region
fh = @(bs) std2(bs.data); % perform this on each block
J = blockproc(IM, bss, fh,'UseParallel',true);
sd=J'; % Column vector
ReportMessage(app,num2str(toc(tstart)));
Is there something I am doing wrong? Do I need to close and reopen the parpool everytime or something.
Also, if I want the fh function to be a functio that has more than 1 argument, how do I do the syntax.
For example, I have a function (that works fine by itself)
Brenner=fmeasure(Image,'BREN',[]); %Focus metric
And I have tried to call it into Blockproc via:
fh = @(bs) fmeasure(bs.data,'BREN',[]);
It doesn't work properly.
1 Comment
"Also, if I want the fh function to be a functio that has more than 1 argument, how do I do the syntax."
If you mean this function:
then it works without error using the syntax you showed:
im = rand(128,4700);
s1 = blockproc(im,[128,20],@(bs)fmeasure(bs.data,'BREN',[]))
Accepted Answer
I doubt that your example is good use of the parallel toolbox. Parallel code can be slower than serial code due to high communication overhead, small task sizes, or data serialization. It may be inefficient when workers send data back to the client frequently or when the overhead of starting a parallel pool outweighs the computation time. And your example has a small task size.
"...analysis of an image (4700 v 128 pixels) that needs breaking up into smaller regions (20x128 pixels)..."
RESHAPE and STD would be simpler and likely faster:
im = rand(128,4700);
tic
jm = reshape(im,128,20,[]);
s0 = std(jm,0,[1,2]);
toc
tic
s1 = blockproc(im,[128,20],@(bs)std2(bs.data));
toc
isequal(s0(:),s1(:))
A FOR loop could be useful too:
for k = 1:size(jm,3)
m = jm(:,:,k);
% do whatever with m
end
... which can be trivially converted to a PARFOR:
34 Comments
"why is that, I thought Blockproc was the ultimate - obviously not.!"
BLOCKPROC splits up data, calls functions repeatedly, works with relatively complex data types, and would have to have a lot of input/error handling. Based on that alone I would assume that it is not very fast. Convenient perhaps, but not fast.
"Could it be the use of std rather than std2?"
I doubt it. It is understanding how data is stored, avoiding splitting data up, avoiding calling functions lots of times, avoiding parallel toolbox for small tasks, using RESHAPE to make data convenient to access, maximising speed by performing multiple operations with single (inbuilt, highly optimised) function calls by using their options, and various other things that one learns. These are the native and natural ways to use MATLAB efficiently.
When a user comes onto this forum as writes "my code is slow and I tried to speed it up using the parallel toolbox" then they have usually misunderstood MATLAB as well as what parallel processing requires and when to apply it.
Jason
9 minuten ago
im = rand(128,4700);
jm = reshape(im,128,20,[]);
s0 = std(jm,0,[1,2]);
size(s0)
sd = s0(:)
hold on
plot(sd)
plot(1/sqrt(12)*ones(size(sd)))
hold off
Actually, im really sorry but Im still running into an error using the reshape approach
My actual image is 4730 pixels wide (not 4700) and 128 pixels high.
I want each sub region to be 20 pixels wide and 128 pixels high like in the picture below
I tried this:
[sy,sx]=size(IM);
nRows=floor(sy/128);
nCols=floor(sx/20);
jm = reshape(IM,nRows,nCols,[]); % reshape = rows, then cols
Error using reshape
Product of known dimensions, 236, not divisible into total number of elements, 605568.
You must take nRows and nCols such that their product gives 4730.
factor(4730)
Thus nRows = 2, nCols = 5*11*43 or nRows = 2*5, nCols = 11*43 or nRows = 5*11, nCols = 2*43 or ... is allowed.
Or you must cut some part of the image.
Jason
ongeveer 2 uur ago
The error message states:
"Product of known dimensions, 236, not divisible into total number of elements, 605568."
Note that the function RESHAPE simply changes the dimension sizes of an array, without adding or removing any data.
Now consider your image, whch contains NP many pixels (where NP does not match your comment):
im = rand(128,4730);
sz = size(im);
np = prod(sz)
Now lets look at what you attempted with those NP pixels:
[sy,sx] = size(im);
nRows = floor(sy/128)
nCols = floor(sx/20)
So you attempted to reshape NP elements into an array with 1 row, 236 columns, and some unknown number of pages. But is NP divisble by 1 and 236 ? Well, 1 is trivial, but 236 lets check if there is any remainder:
mod(np,nCols)
Nope, the number of rows and columns that you attempted to reshape IM into will not work because they cannot evenly divide the number of pixels of IM. That is the cause of the error: RESHAPE cannot create the size array that you requested, because it would require either inventing more data or removing it, which RESHAPE does not do.
Also note that your concept of how to specify the number of rows and columns is unlikely to help you, as should be clear when you perform some basic sanity check, e.g. by asking yourself "why am I reshaping all of the pixels into one row?". You should be counting pixels.
"If i wanted to calculate the Brenner metric instead of the std, how would I do that. It seens std allows for "paging" but its not obvious how the fmeasure function allows for it."
I would skip FMEASURE entirely. If speed is important to you then I would start with a simple FOR loop, which is likely to give reasonable speed, give you a good baseline, and be easy to implement:
im = rand(128,4700);
tic
rn = 128;
cn = 20;
jm = reshape(im,rn,cn,[]);
pn = size(jm,3);
b0 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
a = m(3:end,:)-m(1:end-2,:);
b = m(:,3:end)-m(:,1:end-2);
a(rn,1) = 0;
b(1,cn) = 0;
b0(k) = mean2(max(a,b).^2);
end
toc
b0
"where does this come from?"
It is the expected standard deviation of a continuous uniform distribution U(a,b) between a=0 and b=1. The variance of the uniform distribution is 
which for those a & b gives 
. The standard deviation is therefore just the square root of that.
which for those a & b gives . The standard deviation is therefore just the square root of that.1) because a) inline code tends to be slightly faster than calling a function (as explained here), b) the code for Brenner's metric is fairly simple (just a few lines), c) you can avoid third-party code dependency, d) you may be able to improve/tailor the implementation, e) you learn something about writing code.... and maybe some other reasons. One code-writing rule-of-thumb is do not reinvent the wheel, but this is balanced aginst other rules-of-thumb, such that if you have different priorities (e.g. speed) then tighter coupling actually becomes more useful. In this case there really is nothing stopping you from replacing the entire content of the FOR-loop with an FMEASURE call (which then gets called multiple times).
2) Which part of the loop? The loop preallocates an output vector b0 with NaNs:
then loops over each page of the RESHAPEd image data, storing it in the temporary variable m. Each m is a matrix with size 128x20, which AFAICT is the subimage that you want to process. To calculate Brenner's metric I just looked at the code in FMEASURE, tweaked it, and checked it with an AI tool (btw, you can also look inside FMEASURE). The end metric is Brenner's metric. The metric value is stored in the preallocated vector, there is a vector of them because you have multiple subimages.
But you should not take my word for it: check everything yourself! You can check the size of m yourself using SIZE. You can perform a sanity check on the content of a subimage by extracting the last subimage yourself:
w = im(:,end-19:end);
isequal(m,w)
You can perform a sanity check on the metric by comparing the last metric with the output from FMEASURE:
b1 = fmeasure(m,'BREN',[])
isequal(b0(end),b1)
You can learn more about FOR loops here:
You can practice FOR loops here:
Jason
ongeveer 23 uur ago
Using this function here:
and running on this forum I get:
im = rand(128,4700);
% reshape
rn = 128;
cn = 20;
jm = reshape(im,rn,cn,[]);
pn = size(jm,3);
b0 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
a = m(3:end,:)-m(1:end-2,:);
b = m(:,3:end)-m(:,1:end-2);
a(rn,1) = 0;
b(1,cn) = 0;
b0(k) = mean2(max(a,b).^2);
end
% fmeasure
b1 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
b1(k) = fmeasure(m,'BREN',[]);
end
% comparison:
b0
b1
isequal(b0,b1)
Jason
ongeveer 22 uur ago
It occured to me that you might be providing integer class image data for your actual data, in which case convert it to DOUBLE before running the FOR loop, for example:
jm = reshape(double(im),rn,cn,[]);
I just did some reading about Brenner's metric, and it seems that the implementation in FMEASURE is flawed due to its MAX after squaring. It should square the values before the MAX (or adding them, or whatever).
Consider this carefully-picked example:
im = [10,0,10; 0,0,0; 5,0,5]
fmeasure(im,'BREN',[]) % ZERO !!!!
This occurs because the MAX silently ignores large negative values:
% DH(:,1) = Image(:,3) - Image(:,1)
DH(:,1) = [10-10, 0-0, 5-5] = [0, 0, 0] % First col only; rest stay 0
% → DH is all zeros
% DV(1,:) = Image(3,:) - Image(1,:)
DV(1,:) = [5-10, 0-0, 5-10] = [-5, 0, -5] % First row only; rest stay 0
% → DV has real, significant gradients — but NEGATIVE
FM = max(DH, DV)
= max(0, -5) = 0 % ← the bug! -5 gets silently ignored
FM = mean2(0.^2) = 0
To avoid this square the values before the MAX:
a = im(3:end,:)-im(1:end-2,:);
b = im(:,3:end)-im(:,1:end-2);
a(3,1) = 0;
b(1,3) = 0;
mean2(max(a,b).^2) % wrong!
mean2(max(a.^2,b.^2)) % right!
Some other divergences of FMEASURE from Santos et al:
- it does not implement the gradient threshold v,
- it uses MEAN2 in place of summation.
So I guess that adds another point to the list: f) because random code from the internet is not always reliable: you always need to check it thoroughly.
"This has foxed me!"
Step-by-step debugging: print critical values from each approach, work backwards until you pinpoint where they diverge.
Torsten
ongeveer 22 uur ago
It occured to me that you might be providing integer class image data
im = rand(128,4700);
?
Nick
ongeveer 16 uur ago
Your implementation looks different from Brenner's row version:
function fm = brenner_measure(img)
% Brenner's formula: sum of squared differences 2 pixels apart
diffs = (img(1:end-2, :) - img(3:end, :)).^2;
fm = mean(diffs(:));
end
When you come across negative row and/or column differences your method can break down.
You may try to replace
b0(k) = mean2(max(a,b).^2);
with
b0(k) = mean2(max([a,b,0]).^2);
% or if you want to get upper bound
b0(k) = mean2(max([a,-a,b,-b]).^2);
and compare results
"(Are you saying the formula that Stephen kindly gave me is wrong)?"
I did not give you any formula! Exactly as I wrote in this comment here, I copied the algorithm from FMEASURE. And FMEASURE quotes Santos et al., which I also downloaded. You can find it here:
You selected FMEASURE (which uses Santos et al's algorithm, although perhaps with some bugs), which my code exactly replicates (as demonstrated here).
"Your implementation looks different from Brenner's row version"
Of course it looks different, because it is different! Santos et al explicitly stated their text: "The first three functions (
, 
, 
) were slightly modified to consider both the horizontal and the vertical gradients". Brenner only considered the horizontal gradients, not the horiztonal gradients. But Santos et al stated that they modified this algorithm!
, , ) were slightly modified to consider both the horizontal and the vertical gradients". Brenner only considered the horizontal gradients, not the horiztonal gradients. But Santos et al stated that they modified this algorithm!If you do not even know what algorithm you are using, how do you expect to verify its output?
Have you read Santos et al's paper yet?
"You may try to replace b0(k) = mean2(max(a,b).^2); with b0(k) = mean2(max([a,b,0]).^2);"
No, that won't work because it attempts to concatenate arrays with incompatible sizes. It also attempts to implement something quite different from what Santos et al specified.
" I spent all day trying to track down why Im not getting the same values as you but can't see anything wrong."
You will not find anything by looking all day at some code hoping to see something "wrong". That is not how debugging works.
One obvious start is to compare intermediate values. Have you done that? It should take about two minutes.
Jason
ongeveer 14 uur ago
Edited: Walter Roberson
ongeveer 9 uur ago
Nick
ongeveer 13 uur ago
sorry for sloppy comment — only idea, the code was not correct, just illustration…
a few details follow
for k = 1:pn
m = jm(:,:,k);
a = abs(m(3:end,2:end-1)-m(1:end-2,2:end-1)); % get absolute values
b = abs(m(2:end-1,3:end)-m(2:end-1,1:end-2)); % the same for columns
% a(rn,1) = 0;
% b(1,cn) = 0;
b0(k) = mean2(max(a,b).^2); % now it is correct, but not Brenner
% b0(k) = mean2(max([a,b,0]).^2); — my fault, was careless :
% were not scalar type
end
The original Brenner resembles function brenner_measure() in my previous post
Note that Nicks code discards data without warning. It also has the wrong denominator inside MEAN2.
Consider an image which has gradient only in the first row:
format long G
rn = 3;
cn = 3;
m = [0,4,8; 4,4,4; 4,4,4]
0) Inline code, fixed FMEASURE bug:
a = m(3:end,:)-m(1:end-2,:);
b = m(:,3:end)-m(:,1:end-2);
a(rn,1) = 0;
b(1,cn) = 0;
fm0 = mean2(max(a.^2,b.^2))
1) inline code, same bug as FMEASURE:
fm1 = mean2(max(a,b).^2)
2) FMEASURE
fm2 = fmeasure(m,'BREN',[])
3) Not like anything else:
a = abs(m(3:end,2:end-1)-m(1:end-2,2:end-1)); % get absolute values
b = abs(m(2:end-1,3:end)-m(2:end-1,1:end-2)); % the same for columns
fm3 = mean2(max(a,b).^2)
So even though there is a valid gradient Nick's code does not detect it. The incorrect denominator is slightly more subtle, but is left as an exercise.
Jason
ongeveer 11 uur ago
"Now fmeasure way - need to define fmeasure (with BUG Corrected)"
No, that is not the "fmeasure way", it has at least two bugs:
1) using the entire image to preallocate instead of using zeros, so you are performing all following operations on different values than those which FMEASURE uses.
2) using DIFF which calculates the difference between adjacent elements, not ones separated by some arbitrary step size. As the DIFF documentation clearly states, the 2nd input specifies how many times that this is applied recursively, not some arbitrary step-size. DIFF does not have a step-size option.
After fixing those bugs the code produces exactly the same results as my FOR-loop and the same as the bug-fixed FMEASURE :
im = rand(128,4700);
% reshape
rn = 128;
cn = 20;
jm = reshape(im,rn,cn,[]);
pn = size(jm,3);
% inline
b0 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
a = m(3:end,:)-m(1:end-2,:);
b = m(:,3:end)-m(:,1:end-2);
a(rn,1) = 0;
b(1,cn) = 0;
b0(k) = mean2(max(a.^2,b.^2));
end
% function (fixed MAX bug):
b1 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
b1(k) = fmeasure1(m,'BREN',[]);
end
% code FIXED by me:
b2 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
b2(k) = fmeasureFixed(m);
end
% code BROKEN:
b3 = nan(pn,1);
for k = 1:pn
m = jm(:,:,k);
b3(k) = fmeasureBroken(m);
end
% Compare
b0 % inline
b1 % FMEASURE
b2 % code FIXED by me
b3 % code BROKEN
function FM = fmeasureFixed(Image)
[M,N] = size(Image);
DH = 0*Image; % !!! fixed bug !!!
DV = 0*Image; % !!! fixed bug !!!
DV(1:M-2,:) = Image(3:end,:)-Image(1:end-2,:); % !!! fixed bug !!!
DH(:,1:N-2) = Image(:,3:end)-Image(:,1:end-2); % !!! fixed bug !!!
FM = mean2(max(DH.^2,DV.^2));
end
function FM=fmeasureBroken(Image)
[M,N] = size(Image);
DH = Image; % buggy
DV = Image; % buggy
DH(1:M-2,:) = diff(Image,2,1); % buggy
DV(:,1:N-2) = diff(Image,2,2); % buggy
FM = mean2(max(DH.^2,DV.^2));
end
Also note that your comment "Bug where -ve values are ignored." is incorrect. I never wrote that all negative values are ignored.
Jason
ongeveer 4 uur ago
Hi Stephen, you've done it!! You've worked out why I wasn't matching what you had. Its because the version of fmeasure I was using was pre 2016 .... where DH=image and the diff function was used. You must be using a latest version
One last query (I understand if you've had enough now - and I really appreciate the massive effort you have put into this).
How can "some" -ve values be ignored and not all?
"How can "some" -ve values be ignored and not all?"
max(-3,-9)
Nick
21 minuten ago
Resume:
@Stephen23 analyzed codes for "Brenner's" variant of fmeasure() function:
0) buggy code: wrong results as a rule; better forget it
2) no code to analyze
3) corrected by skipping boundary gradients, looks like not what expected
To use a more robust version of the function with boundaries taken into account, try:
m = [0,4,8; 4,4,4; 4,4,4; 4,0,4]'
a = m(3:end,:)-m(1:end-2,:);
b = m(:,3:end)-m(:,1:end-2);
fm4 = mean2([a(:).^2;b(:).^2])
P.S. The code only demonstrates an approach to get "mean gradient" inside domain (without the boundary). That is why no warnings and error handlers assumed. BTW denominator is correct.
More Answers (0)
Categories
Find more on Tables in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
