我正在尝试实现本文描述的算法:
临时光谱带中生物斑点图像的分解
这是算法的说明:
我们N以采样频率记录了一系列连续的斑点图像fs。这样,可以观察像素如何在图像中演变N。可以将这种演变视为时间序列,并可以通过以下方式进行处理:与每个像素的演变相对应的每个信号都用作一组滤波器的输入。预先将强度值除以它们的时间平均值,以最小化对象的反射率或照明的局部差异。可以充分分析的最大频率取决于采样定理和采样频率的一半fs。后者由CCD相机,图像尺寸和图像采集卡设置。滤波器组如图1所示。 在我们的案例中,使用了十个5°阶Butterworth滤波器,但是该数量可以根据所需的区分度而变化。该银行是使用MATLAB软件在计算机上实现的。我们选择Butter- worth滤波器是因为它不仅简单,而且最大程度地平坦。可以使用其他滤波器,无限冲激响应或有限冲激响应。 通过该组滤波器,获得了每个临时像素演变的每个滤波器的十个对应信号作为输出。然后计算每个信号中的平均能量Eb: 其中,pb(n)是第n个要过滤的图像中经过滤波的像素的强度b除以其平均值,并且N是图像总数。这样,En获得了每个像素的能量值,每个下摆都属于图1中的一个频带。 利用这些值,可以构建活动对象的十个图像,每个图像都显示在某个频带中存在多少随时间变化的散斑能量。错误地将颜色分配给结果中的灰度级将有助于区分。
我们N以采样频率记录了一系列连续的斑点图像fs。这样,可以观察像素如何在图像中演变N。可以将这种演变视为时间序列,并可以通过以下方式进行处理:与每个像素的演变相对应的每个信号都用作一组滤波器的输入。预先将强度值除以它们的时间平均值,以最小化对象的反射率或照明的局部差异。可以充分分析的最大频率取决于采样定理和采样频率的一半fs。后者由CCD相机,图像尺寸和图像采集卡设置。滤波器组如图1所示。
N
fs
在我们的案例中,使用了十个5°阶Butterworth滤波器,但是该数量可以根据所需的区分度而变化。该银行是使用MATLAB软件在计算机上实现的。我们选择Butter- worth滤波器是因为它不仅简单,而且最大程度地平坦。可以使用其他滤波器,无限冲激响应或有限冲激响应。
通过该组滤波器,获得了每个临时像素演变的每个滤波器的十个对应信号作为输出。然后计算每个信号中的平均能量Eb:
其中,pb(n)是第n个要过滤的图像中经过滤波的像素的强度b除以其平均值,并且N是图像总数。这样,En获得了每个像素的能量值,每个下摆都属于图1中的一个频带。
pb(n)
b
En
利用这些值,可以构建活动对象的十个图像,每个图像都显示在某个频带中存在多少随时间变化的散斑能量。错误地将颜色分配给结果中的灰度级将有助于区分。
这是我基于此的MATLAB代码:
for i=1:520 for j=1:368 ts = []; for k=1:600 ts = [ts D{k}(i,j)]; %%% kth image pixel i,j --- ts is time series end ts = double(ts); temp = mean(ts); if (temp==0) for l=1:10 filtImag1{l}(i,j)=0; end continue; end ts = ts-temp; ts = ts/temp; N = 5; % filter order W = [0.0 0.10;0.10 0.20;0.20 0.30;0.30 0.40;0.40 0.50;0.50 0.60 ;0.60 0.70;0.70 0.80 ;0.80 0.90;0.90 1.0]; [B,A]=butter(N,0.10,'low'); ts_f(1,:) = filter(B,A,ts); N1 = 5; for ind = 2:9 Wn = W(ind,:); [B,A] = butter(N1,Wn); ts_f(ind,:) = filter(B,A,ts); end [B,A]=butter(N,0.90,'high'); ts_f(10,:) = filter(B,A,ts); for ind=1:10 %Following Paper Suggestion filtImag1{ind}(i,j) =sum(ts_f(ind,:).^2); end end end for i=1:10 figure,imshow(filtImag1{i}); colorbar end pre_max = max(filtImag1{1}(:)); for i=1:10 new_max = max(filtImag1{i}(:)); if (pre_max<new_max) pre_max=max(filtImag1{i}(:)); end end new_max = pre_max; pre_min = min(filtImag1{1}(:)); for i=1:10 new_min = min(filtImag1{i}(:)); if (pre_min>new_min) pre_min = min(filtImag1{i}(:)); end end new_min = pre_min; %normalize for i=1:10 temp_imag = filtImag1{i}(:,:); x=isnan(temp_imag); temp_imag(x)=0; t_max = max(max(temp_imag)); t_min = min(min(temp_imag)); temp_imag = (double(temp_imag-t_min)).*((double(new_max)-double(new_min))/double(t_max-t_min))+(double(new_min)); %median filter %temp_imag = medfilt2(temp_imag); imag_test2{i}(:,:) = temp_imag; end for i=1:10 figure,imshow(imag_test2{i}); colorbar end for i=1:10 A=imag_test2{i}(:,:); B=A/max(max(A)); B=histeq(A); figure,imshow(B); colorbar imag_test2{i}(:,:)=B; end
但我得到的结果与纸张不同。有谁知道为什么?还是我哪里出问题了?
*通过从@Amro获得帮助并使用他的代码进行 *编辑 ,我得到了以下图像:这是我从72hrs发芽的扁豆中提取的原始图像(400张图像,每秒5帧):
这是10个不同波段的结果图像:
我可以发现几个问题:
当您将信号除以均值时,需要检查它是否不为零。否则结果将是NaN。
NaN
作者(我正在关注本文)使用了一组滤波器,这些滤波器的频带覆盖了直到奈奎斯特频率的整个范围。您正在做一半。您传递给的归一化频率butter应该一直上升到1(对应于fs/2)
butter
1
fs/2
在计算每个滤波信号的能量时,我认为您不应该除以其均值(您之前已经考虑过了)。取而代之的是:E = sum(sig.^2);对每个滤波信号
E = sum(sig.^2);
在最后的后处理步骤中,应将其归一化为[0,1]范围,然后应用中值滤波算法medfilt2。计算看起来不正确,应该是这样的:
medfilt2
img = ( img - min(img(:)) ) ./ ( max(img(:)) - min(img(:)) );
考虑到以上几点,我试图以向量化的方式重写代码。由于您没有发布示例输入图像,因此我无法测试结果是否符合预期。此外,我不确定如何解释最终图像:)
%# 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)
(我将 此处的图片用于上述结果)
进行了一些更改并简化了上面的代码。这将减少内存占用。例如,我使用单元格数组而不是单个多维矩阵来存储结果。这样,我们就不会分配一大块连续内存。我还重用了相同的变量,而不是在每个中间步骤中都引入了新变量。