StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>A couple of issue I can spot:</p> <ul> <li><p>when you divide the signal by its mean, you need to check that it was not zero. Otherwise the result will be <code>NaN</code>.</p></li> <li><p>the authors (I am following <a href="http://wheat.pw.usda.gov/ggpages/BarleyNewsletter/50/ARGBNL50.htm" rel="nofollow noreferrer">this article</a>) used a bank of filters with frequency bands covering the entire range up to the Nyquist frequency. You are doing half of that. The normalized frequencies you pass to <code>butter</code> should go all the way up to <code>1</code> (corresponds to <code>fs/2</code>)</p></li> <li><p>When computing the energy of each filtered signal, I think you should not divide by its mean (you have already accounted for that before). Instead simply do: <code>E = sum(sig.^2);</code> for each of the filtered signals</p></li> <li><p>In the last post-processing step, you should normalize to the range [0,1], and then apply the median filtering algorithm <code>medfilt2</code>. The computation doesn't look right, it should be something like:</p> <pre><code>img = ( img - min(img(:)) ) ./ ( max(img(:)) - min(img(:)) ); </code></pre></li> </ul> <hr> <h3>EDIT:</h3> <p>With the above points in mind, I tried to rewrite the code in a vectorized way. Since you didn't post sample input images, I can't test if the result is as expected... Plus I am not sure how to interpret the final images anyway :)</p> <pre><code>%# read biospeckle images fnames = dir( fullfile('folder','myimages*.jpg') ); fnames = {fnames.name}; N = numel(fnames); %# number of images Fs = 1; %# sampling frequency in Hz sz = [209 278]; %# image sizes T = zeros([sz N],'uint8'); %# store all images for i=1:N T(:,:,i) = imread( fullfile('folder',fnames{i}) ); end %# timeseries corresponding to every pixel T = reshape(T, [prod(sz) N])'; %# columns are the signals T = double(T); %# work with double class %# normalize signals before filtering (avoid division by zero) mn = mean(T,1); T = bsxfun(@rdivide, T, mn+(mn==0)); %# divide by temporal mean %# bank of filters numBanks = 10; order = 5; % butterworth filter order fCutoff = linspace(0, Fs/2, numBanks+1)'; % lower/upper cutoff freqs W = [fCutoff(1:end-1) fCutoff(2:end)] ./ (Fs/2); % normalized frequency bands W(1,1) = W(1,1) + 1e-5; % adjust first freq W(end,end) = W(end,end) - 1e-5; % adjust last freq %# filter signals using the bank of filters Tf = cell(numBanks,1); %# filtered signals using each filter for i=1:numBanks [b,a] = butter(order, W(i,:)); %# bandpass filter Tf{i} = filter(b,a,T); %# apply filter to all signals end clear T %# cleanup unnecessary stuff %# compute average energy in each signal across frequency bands Tf = cellfun(@(x)sum(x.^2,1), Tf, 'Uniform',false); %# normalize each to [0,1], and build corresponding images Tf = cellfun(@(x)reshape((x-min(x))./range(x),sz), Tf, 'Uniform',false); %# show images for i=1:numBanks subplot(4,3,i), imshow(Tf{i}) title( sprintf('%g - %g Hz',W(i,:).*Fs/2) ) end colormap(gray) </code></pre> <p><img src="https://i.stack.imgur.com/bBoFB.png" alt="screenshot"></p> <p><em>(I used the image from <a href="http://lissi.fr/eric-delechelle/research/speckle/start" rel="nofollow noreferrer">here</a> for the above result)</em></p> <h3>EDIT#2</h3> <p>Made some changes and simplified the above code a bit. This shall reduce memory footprint. For example I used cell array instead of a single multidimensional matrix to store the result. That way we don't allocate one big block of contiguous memory. I also reused same variables instead of introducing new ones at each intermediate step...</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload