2

this post is related to my previous question : image processing in matlab as I have uploaded my algorithme there. the think is that I am trying to change the filtering part of the code. in matlab filter.m function can accpet filter(B, A, my pixels evolutions in Time) and it return me the filtered values. but at the moment I have to pass the the whole time series together.

but the problem is that now I want to change the code in a way instead of passing the whole timeseries into the filter, I want to pass one value at a time, but I want filter function treat the value like the nth value not the first valu. I have created a sudo code, as I am injecting one picture into the code, but I dont know how can change the filtering part., any body has any idea??

clear all
j=1;
for i=0:3000
  str = num2str(i);
  str1 = strcat(str,'.mat');
  load(str1);
  D{j}=A(20:200,130:230);
  j=j+1;
end
N=5;
w = [0.00000002 0.05;0.05 0.1;0.1 0.15;0.15 0.20;
     0.20 0.25;0.25 0.30;0.30 0.35;0.35 0.40;
     0.40 0.45;0.45 0.50;0.50 0.55;0.55 0.60;
     0.60 0.65;0.65 0.70;0.70 0.75;0.75 0.80;
     0.80 0.85;0.85 0.90;0.90 0.95;0.95 0.99999999];
for ind=1:20
  wn = w(ind,:);
  [b,a] = butter(N,wn);
  bCoeff(ind,:)=b;
  aCoeff(ind,:)=a;
end

ts=[];
sumLastImages=[];
for k=1:10 %number of images
  for bands=1:20 %number of bands
    for i=1:10 %image h
      for j=1:10 %image w
        pixelValue = D{k}(i,j);          
        % reflectivity elimination        
        % for the current pixel, have the summation of the same position from before 
        % images and create a mean value base on the temporal values
        sumLastImages(i,j)=pixelValue+sumLastImages(i,j);
        meanValue = sumLastImages(i,j)/k;

        if(meanValue==0)            
          filteredimages{bands}(i,j)=0; 
          continue;
        else
          pixel_calculated_meanvalue = pixelValue/meanValue; 
        end                      
        % filter part that need to be changed, and because we are adding 
        % one value then the reutrn of the filter is one too         
        ts_f = filter(bCoeff(bands,:), aCoeff(bands,:), ...
                      pixel_calculated_meanvalue);                         
        filteredimages{bands}(i,j)=ts_f;            
      end     
    end       

    finalImagesSummation{bands}(:,:) = ...
      (filteredimages{bands}(:,:)^2)+finalImagesSummation{bands}(:,:);
    finalImages{bands}(:,:)=finalImagesSummation/k;
  end
end

EDIT I managed to change the code like this, which now I load the fist 200 images, and after that I am able to inject the images one by one into the filter, but now the problem is that I am getting "Out of memory. Type HELP MEMORY for your options." error for big images. here is my code any idea to efficent the code :

%%
cd('D:\MatlabTest\06-06-Lentils');
clear all

%%

N=5;
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];

[bCoeff{1},aCoeff{1}] = butter(N,0.1,'Low');
for ind=2:9
  wn = W(ind,:);
  [b,a] = butter(N,wn);
  bCoeff{ind}=b;
  aCoeff{ind}=a;
end
[bCoeff{10},aCoeff{10}] = butter(N,0.9,'high');

%%
j=1;

D = zeros(200,380,320);

T = 200;
K = 0:399;
l = T+1;
Yout = cell(1,10);
figure;

for i = 1:length(K)-200
  disp(i)

  if i == 1
    for j = 1:T
      str = int2str(K(1)+j-1);
      str1 = strcat(str,'.mat');
      load(str1);
      D(j,1:380,1:320) = A;
    end

  else

    str = int2str(l);
    str1 = strcat(str,'.mat');
    load(str1);

    temp = D(2:200,1:380,1:320) ;
    temp(200,1:380,1:320) = A;
    D = temp;
    clear temp
    l = l +1;

  end


  for p = 1:380
    for q = 1:320
      x = D(:,p,q) - mean(D(:,p,q));

      for k = 1:10
        temp = filter(bCoeff{k},aCoeff{k},x);
        if i == 1
          Yout{k}(p,q) = mean(temp);
        else
          Yout{k}(p,q) = (Yout{k}(p,q)  + mean(temp))/2;
        end
      end
    end
  end

  for k = 1:10
    subplot(5,2,k)
    subimage(Yout{k}*1000,[0 255]);
    color bar
    colormap jet
  end
  pause(0.01);
end

disp('Done Loading...')
Community
  • 1
  • 1
user261002
  • 2,182
  • 10
  • 46
  • 73
  • I think you should start at the basic level: How many images do you want to load in parallel and how large are they (width, height, bits per pixel)? – Mehrwolf Jul 25 '12 at 09:28
  • @Mehrwolf the stack of images and the size always vary, but the camera is a monochromic camera 8bits per pixel – user261002 Jul 25 '12 at 09:51
  • But you have a rough estimate for the amount of images and their size, right? – Mehrwolf Jul 25 '12 at 11:08
  • Another thing: If the image has only 8 bits per pixel, store it in a uint8 array (type `help uint8`). Matlab uses double arrays per default, which is 8 times as large. – Mehrwolf Jul 25 '12 at 11:10

3 Answers3

2

Overview

If all you want to have is an IIR filter, which you can feed incrementally, i.e. where you do not have to supply the full vector at once, you can simply implement your own function and use either persistent variables.

Simple approach

The code snippet below defines a function myFilter, which makes use of persistent variables, which you can control with the following commands:

  • init: set up the IIR coefficients
  • getA: returns the A coefficients, i.e. the outputweights
  • getB: returns the B coefficients, i.e. the input weights
  • getX: returns the stored input data x[0], x[1], ... x[M]
  • getY: returns the output data y[0], y[1], ... y[N]
  • getCurrentY: returns the last output data y[N]

Here is the function:

function result = myFilter(varargin)
% myFilter   A simple IIR filter which can be fed incrementally.
%
% The filter is controlled with the following commands:
%   myFilter('init', B, A)
%      Initializes the coefficients B and A. B are the weights for the
%      input and A for the output.
%   myFilter('reset')
%      Resets the input and output buffers to zero.
%   A = myFilter('getA')
%   B = myFilter('getB')
%      Returns the filter coefficients A and B, respectively.
%   x = myFilter('getX')
%   y = myFilter('getY')
%      Returns the buffered input and output vectors.
%   y = myFilter('getCurrentY')
%      Returns the current output value.
%   myFilter(x)
%      Adds the new value x as input to the filter and updates the 
%      output.

    persistent Bcoeff
    persistent Acoeff
    persistent x
    persistent y

    if ischar(varargin{1})
        % This is a command.
        switch varargin{1}
            case 'init'
                Bcoeff = varargin{2};
                Acoeff = varargin{3};
                Bcoeff = Bcoeff / Acoeff(1);
                Acoeff = Acoeff / Acoeff(1);
                x = zeros(size(Bcoeff));
                y = zeros(1, length(Acoeff) - 1);
                return

            case 'reset'
                x = zeros(size(Bcoeff));
                y = zeros(1, length(Acoeff) - 1);
                return

            case 'getA'
                result = Acoeff;
                return

            case 'getB'
                result = Bcoeff;
                return

            case 'getX'
                result = x;
                return

            case 'getY'
                result = y;
                return

            case 'getCurrentY'
                result = y(1);
                return

            otherwise
                error('Unknown command');
        end
    end

    % A new value has to be filtered.
    xnew = varargin{1};
    x = [xnew, x(1:end-1)];
    ynew = sum(x .* Bcoeff) - sum(y .* Acoeff(2:end));
    y = [ynew, y(1:end-1)];
end

And a usage example:

% Define the filter coefficients. Bcoeff acts on the input, Acoeff on
% the output.
Bcoeff = [4, 5];
Acoeff = [1, 2, 3];
% Initialize the filter.
myFilter('init', Bcoeff, Acoeff);

% Add a value to be filtered.
myFilter(10)
myFilter('getCurrentY')

ans =
    40

% Add another one.
myFilter(20)
myFilter('getCurrentY')

ans =
    50

% And a third one.
myFilter(30)
myFilter('getCurrentY')

ans =
     0

% Compare with the Matlab filter function.
filter(Bcoeff, Acoeff, [10 20 30])

ans =
    40    50     0

The drawback of this approach is that it is only possible to have one active filter simultaneously. This is problematic e.g. in your question, where you have different filters which are updated in an alternating fashion.

Advanced approach

In order to operate multiple filters simultatenously, you need some way to identify the filter. The solution I present here works with handles. A handle is simple an integer. To be more precise, it is actually an index into a persistent array, which itself holds the filter state, i.e. the filter coefficients and the buffers for the input and the output.

The calling syntax is a bit more complicated, because you have to pass a handle. However, it is possible that multiple filters are active simultaneously.

function result = myFilterH(varargin)
% myFilterH   A simple IIR filter which can be fed incrementally.
%             Operates on a filter handle.
%
% The filter is controlled with the following commands:
%   handle = myFilterH('create')
%      Creates a new filter handle.
%   myFilterH(handle, 'init', B, A)
%      Initializes the coefficients B and A. B are the weights for the
%      input and A for the output. handle identifies the filter.
%   myFilterH(handle, 'reset')
%      Resets the input and output buffers to zero.
%   A = myFilterH(handle, 'getA')
%   B = myFilterH(handle 'getB')
%      Returns the filter coefficients A and B, respectively.
%   x = myFilterH(handle, 'getX')
%   y = myFilterH(handle, 'getY')
%      Returns the buffered input and output vectors.
%   y = myFilterH(handle, 'getCurrentY')
%      Returns the current output value.
%   myFilterH(handle, x)
%      Adds the new value x as input to the filter and updates the 
%      output.

    persistent handles

    if ischar(varargin{1})
        if strcmp(varargin{1}, 'create')
            if isempty(handles)
                handles = struct('x', [], 'y', [], 'A', [], 'B', []);
                result = 1;
            else
                result = length(handles) + 1;
                handles(result).x = [];
            end
            return
        else
            error('First argument must be a filter handle or ''create''');
        end
    end

    % The first input should be the handles(hIdx).
    hIdx = varargin{1};
    if hIdx < 0 || hIdx > length(handles)
        error('Invalid filter handle')
    end

    if ischar(varargin{2})
        % This is a command.
        switch varargin{2}
            case 'init'
                handles(hIdx).B = varargin{3};
                handles(hIdx).A = varargin{4};
                handles(hIdx).B = handles(hIdx).B / handles(hIdx).A(1);
                handles(hIdx).A = handles(hIdx).A / handles(hIdx).A(1);
                handles(hIdx).x = zeros(size(handles(hIdx).B));
                handles(hIdx).y = zeros(1, length(handles(hIdx).A) - 1);
                return

            case 'reset'
                handles(hIdx).x = zeros(size(handles(hIdx).B));
                handles(hIdx).y = zeros(1, length(handles(hIdx).A) - 1);
                return

            case 'getA'
                result = handles(hIdx).A;
                return

            case 'getB'
                result = handles(hIdx).B;
                return

            case 'getX'
                result = handles(hIdx).x;
                return

            case 'getY'
                result = handles(hIdx).y;
                return

            case 'getCurrentY'
                result = handles(hIdx).y(1);
                return

            otherwise
                error('Unknown command');
        end
    end

    % A new value has to be filtered.
    xnew = varargin{2};
    handles(hIdx).x = [xnew, handles(hIdx).x(1:end-1)];
    ynew = sum(handles(hIdx).x .* handles(hIdx).B) ...
               - sum(handles(hIdx).y .* handles(hIdx).A(2:end));
    handles(hIdx).y = [ynew, handles(hIdx).y(1:end-1)];
end

And the example:

% Define the filter coefficients.
Bcoeff = [4, 5];
Acoeff = [1, 2, 3];
% Create new filter handles.
fh1 = myFilterH('create');
fh2 = myFilterH('create');
% Initialize the filter handle. Note that reversing Acoeff and Bcoeff creates
% two totally different filters.
myFilterH(fh1, 'init', Bcoeff, Acoeff);
myFilterH(fh2, 'init', Acoeff, Bcoeff);

% Add a value to be filtered.
myFilterH(fh1, 10);
myFilterH(fh2, 10);
[myFilterH(fh1, 'getCurrentY'), myFilterH(fh2, 'getCurrentY')]

ans =
    40.0000    2.5000

% Add another one.
myFilterH(fh1, 20);
myFilterH(fh2, 20);
[myFilterH(fh1, 'getCurrentY'), myFilterH(fh2, 'getCurrentY')]

ans =
    50.0000    6.8750

% And a third one.
myFilterH(fh1, 30);
myFilterH(fh2, 30);
[myFilterH(fh1, 'getCurrentY'), myFilterH(fh2, 'getCurrentY')]

ans =
     0   16.4063

% Compare with the Matlab filter function.
filter(Bcoeff, Acoeff, [10 20 30])

ans =
    40    50     0

filter(Acoeff, Bcoeff, [10 20 30])
ans =
    2.5000    6.8750   16.4063

Using it for your example

To use the advanced approach in your example, you could first create an array of filters:

fh = [];
for ind = 1:20
  wn = w(ind,:);
  [b,a] = butter(N,wn);
  fh(ind) = myFilterH('create');
  myFilterH(fh(ind), 'init', b, a);
end

Later on in the filter part simply add your value to all of the filters. However, for this you also need to reverse the loops because right now you would need one filter per band per pixel. If you loop over pixels first and then over bands, you can reuse the filters:

for height = 1:10
  for width = 1:10 
    for bands = 1:20
      myFilterH(fh(bands), 'reset');
      for k = 1:10
        [...]
        ts_f = myFilterH(fh(bands), ...
                         pixel_calculated_meanvalue);                         
        filteredimages{bands}(i,j) = myFilterH(fh(bands), 'getCurrentY');
      end
    end
  end
end

Hope that helps!

Mehrwolf
  • 8,208
  • 2
  • 26
  • 38
  • a third option: nested functions and closures to maintain state between calls instead of using persistent/global variables. – Amro Jul 25 '12 at 14:10
  • @Amro: How could such a thing look like? Would I then need to embed the nested functions and the closures in the code which calls the filter? – Mehrwolf Jul 26 '12 at 09:54
  • @Amro , @Mehrwolf there is actually no need to all of this. `filter` already has this built-in functionality. – Eitan T Jul 26 '12 at 21:12
  • @EitanT: I admit I haven't read this question in detail, I was just looking at it from the point of view of variables scope and maintaining state between calls.. Your solution definitely makes more sense than what I suggested +1 :) – Amro Jul 26 '12 at 21:21
2

No need to rewrite the filter function, there is a simple solution!

If you want to feed filter with one sample at a time, you need to pass the state parameters as well so that each input sample is processed depending on its predecessor. After filtering, the new state is actually returned as a second parameter, so that most of the work is already done by MATLAB for you. This is good news!

For the sake of readability, allow me to temporarily replace your variable names with simple letters:

a = aCoeff(bands, :);    
b = bCoeff(bands, :);
x = pixel_calculated_meanvalue;

ts_f is represented by y.

And so, this:

y = filter(b, a, x);

is actually equivalent to this:

N = numel(x);
y = zeros(N, 1);                             %# Preallocate memory for output
z = zeros(max(length(a), length(b)) - 1, 1); %# This is the initial filter state
for i = 1:N
    [y(i), z] = filter(b, a, x(i), z);
end

You can check for yourself that the result is the same!

For your example, the code would be:

N = numel(pixel_calculated_meanvalue);
ts_f = zeros(N, 1);
z = zeros(max(length(aCoeff(bands, :)), length(bCoeff(bands, :))) - 1, 1); 
for i = 1:N
    [ts_f(i), z] = filter(bCoeff(bands, :), aCoeff(bands, :), ...
        pixel_calculated_meanvalue(i), z);
end

With this method you can process one input sample at a time, just make sure you store the last filter state after every filter call. If you plan on using multiple filters, you'll have to store a state vector per filter!

Eitan T
  • 32,660
  • 14
  • 72
  • 109
0

If I understand the question correctly, the crux is "How do I return multiple things from a Matlab function?"

You can return multiple things like this:

function [a,  b, np, x, y] = filter(ord, a,  b, np, x, y)
    %code of function here
    %make some changes to a, b, np, x, and y
end

If you want to call the filter from another function and catch its return values, you can do this:

function [] = main()
    %do some stuff, presumably generate initial values for ord, a,  b, np, x, y
    [new_a,  new_b, new_np, new_x, new_y] = filter(ord, a,  b, np, x, y)  
end

In short, if you do function [x, y] = myfunc(a, b), then myfunc returns both x and y.

solvingPuzzles
  • 8,541
  • 16
  • 69
  • 112
  • not at all, I want to change the filter implementation in matlab,but I dont know how and I have a c++ code that I dont know if it is correct and can I replace it with filter function in matlab to desing the butterworth filter – user261002 Jul 20 '12 at 16:54