Weird Behaviour of Parallel Processing when Using BlockProc
Mostrar comentarios más antiguos
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 comentario
"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',[]))
Respuesta aceptada
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:
29 comentarios
"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
hace 9 minutos
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
hace alrededor de 1 hora
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
hace alrededor de 10 horas
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
hace alrededor de 9 horas
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
hace alrededor de 8 horas
It occured to me that you might be providing integer class image data
im = rand(128,4700);
?
Nick
hace alrededor de 3 horas
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
hace 29 minutos
Nick
hace alrededor de 2 horas
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
hace 20 minutos
Más respuestas (0)
Categorías
Más información sobre Blocked Images en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
