% *************************************************************************

% ECE431F Digital Signal Processing

%

% Lab 4: Practical application of DSP to a real world problem:

% Processing digital signals from an infrared sensor operated faucet.

%

% Matlab m. file scripts, you can rename it to lab4.m if you like.

%

% Answer the 20 questions in your report.

% **************************************************************************

% ***************************** %

% Introduction and Instructions %

% ***************************** %

% You are required to use DSP to figure out how to control an

% infrared sensor faucet. When hands are placed under the tap

% the controller should turn the water on and when the hands

% are removed the controller should turn the water off.

% There are a number of ways to do this. The method we will use is

% based on image processing using masks and a simple probability

% model of the water and hands. To do this we will need to compute

% means, standard deviations, and masks of useful object in our

% image space.

% The lab is broken up into 4 parts.

% A. Computing Means, Variances, and Masks.

% B. Using an estimate of the total image power

% C. Computing the power intelligently using masks

% D. Using Probability to make more intelligent decisions

% There are four graphs to produce and twenty questions to answer.

% These are not essay questions, they are short answer and should

% be no more than 3 to 5 sentences long (except for Q18 which we'd

% like you to elaborate on). You should be able to do

% this lab in a couple of hours, the code has been largely written

% for you, all you need to do is plot out some of the results and

% understand what is going on.

% Don't be frightened by all the for-loops, much of the code is

% repetitive.

% If you find that you are running out of RAM or your machine is too

% slow you can adjust the range of the for loops within reason.

% Good luck and happy learning!

% ***************************** %

% Preliminary Image Information %

% ***************************** %

% Clear all previous data

clear all;

dispon=0;%display processing information or not

% Image information

M = 240; % Image height in pixels

N = 320; % Image width in pixels

% The sample images are stored in different directories and the directories

% in which these directories are stored varies between users. Please adjust

% the following few lines of code to match your storage locations

rootdir = './';

filedirsink = 'noise0/';

filedirwater = 'noise1/';

filedirhands = 'signal0/';

filedirhandswater = 'signal1/';

filetype = 'jpg';

disp('Part A:');

% *********************************** %

% %

% Part A: Means, Variances, and Masks %

% %

% *********************************** %

% Compute the mean of the Handwash Sink Faucet signal without hands and water

meanA = zeros(M,N);

L = 58; % number of frames in dataset

for i=0:L-1

filename = strcat('a',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirsink,filename),'jpeg'));

if dispon

msg=sprintf('a mean %d\r',i);

disp(msg);

end

meanA = meanA + testimage/L;

end%for

%%% Q1. What does the mean represent?

%%% Q2. What is the advantage of computing the mean of a large number of

%%% realizations of the same signal over just taking one of these signals?

% Compute the standard deviation about the mean of the sink

varA = zeros(M,N);

for i=0:L-1

filename = strcat('a',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirsink,filename),'jpeg'));

if dispon

msg=sprintf('a var %d\r',i);

disp(msg);

end

varA = varA + (double(testimage) - meanA).^2/(L-1);

end%for

stdA = sqrt(varA);

%%% Q3. Is the image of the sink stochastic in nature? Explain.

%%% Q4. What does the standard deviation represent for the sink image?

%%% Q5. Approximately how many standard deviations from the mean do we have

%%% to be to contain 95% of our samples?

%%% Q6. Approximately what percentage of the samples are within 3 standard

%%% deviations of the mean?

% Compute mean when the water is flowing:

meanB = zeros(M,N);

maskB = zeros(M,N);

L=155; % number of measurements in the dataset

for i=0:L-1,

filename = strcat('b',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirwater,filename),'jpeg'));

if dispon

msg=sprintf('b mean %d\r',i);

disp(msg);

end

meanB = meanB + testimage/L;

maskB = maskB + (testimage - meanA)/L;

end%for

threshold = -30;

maskB = (maskB < threshold);

%%% Q7. In this context, what is a mask?

%%% Is it a filter?

%%% What is the minimum and maximum value of maskB?

%%% Q8. Why is it important to subtract meanA from the testimages in the

%%% determination of the mask?

%%% Q9. Plot the image of maskB for 3 different thresholds

%%% (eg.-50,-30, and -10). Is -30 a good threshold? Explain.

% Compute the standard deviation when the water is flowing in the sink

varB = zeros(M,N);

for i=0:154,

filename = strcat('b',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirwater,filename),'jpeg'));

if dispon

msg=sprintf('b var %d\r',i);

disp(msg);

end

varB = varB + (testimage - meanB).^2/154;

end%for

stdB = sqrt(varB);

%%% Q10. Is the image of the water alone stochastic in nature? Explain.

%%% Q11. What does the standard deviation represent for the water image?

% Compute the mean of the hands

meanC = zeros(M,N);

maskC = zeros(M,N);

for i=33:41,

filename = strcat('d',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhands,filename),'jpeg'));

if dispon

msg=sprintf('c mean %d\r',i);

disp(msg);

end

meanC = meanC + testimage/9;

maskC = maskC + (testimage - meanA)/9;

end%for

maskC = (maskC < threshold);

% Compute the standard deviation of the hands

varC = zeros(M,N);

for i=33:41,

filename = strcat('d',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhands,filename),'jpeg'));

if dispon

msg=sprintf('c var %d\r',i);

disp(msg);

end

varC = varC + (testimage - meanC).^2/8;

end%for

stdC = sqrt(varC);

% The final image masks

maskHands = maskC - maskB.*maskC;

maskWater = maskB;

%%% Plot: Figure 1. Use imagesc to plot the following grayscale plots:

%%% [1,1] meanA [1,2] stdA [1,3] maskHands

%%% [2,1] meanB [2,2] stdB [2,3] maskWater

%%% [3,1] meanC [3,2] stdC [2,3] maskC

colormap(gray); % colormaps usually come before imagesc

subplot(3,3,1),imagesc(meanA);colormap(gray);title('1.1. The Sink')

subplot(3,3,2),imagesc(stdA);colormap(gray);title('1.2. Standard Deviation of Sink')

subplot(3,3,3),imagesc(maskHands);colormap(gray);title('1.3. Mask of Hands Only')

subplot(3,3,4),imagesc(meanB);colormap(gray);title('1.4. The Water')

subplot(3,3,5),imagesc(stdB);colormap(gray);title('1.5. Standard Deviation of Water')

subplot(3,3,6),imagesc(maskWater);colormap(gray);title('1.6. Mask of Water Only')

subplot(3,3,7),imagesc(meanC);colormap(gray);title('1.7. The Hands')

subplot(3,3,8),imagesc(stdC);colormap(gray);title('1.8. Standard Deviation of Hands')

subplot(3,3,9),imagesc(maskC);colormap(gray);title('1.9. Mask of Water and Hands')

%%% Q12. What does maskHands represent?

disp('part B')

% ************************************************** %

% %

% Part B: Using L1 norm to control the faucet %

% %

% ************************************************** %

% In this section we use our first intuitive guess of how the faucet control

% should operate based on the images. For this lab we define the L1 norm as

% the cumulative sum of the absolute values of the pixels of an image.

% That is,

%

% sum2(|image|) ~ sum(sum(|image|))

%

% There are other norms (e.g. power which is L2), but this simplistic one

% will suffice for this lab. Whenever the hands or the water appear above

% the sink the image appears to gain L1 norm. So consider the design of a

% simple control, using the power level to determine when to turn the water

% on and when to turn it off.

L=58; % Length of image sequence

Pa = NaN*ones(L,1); % initialize the column vector of L1 'probabilities'

% this is not necessary but it's more efficient

for i=0:L-1,

filename = strcat('a',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirsink,filename),'jpeg'));

Pa(i+1) = sum(sum(abs(testimage - meanA)));

end%for

L=155; % Length of image sequence

Pb = NaN*ones(L,1); % initialize the column vector of L1 'probabilities'

for i=0:L-1,

filename = strcat('b',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirwater,filename),'jpeg'));

Pb(i+1) = sum(sum(abs(testimage - meanA)));

end%for

P =[];

for i=24:41,

filename = strcat('d',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhands,filename),'jpeg'));

P = [P; sum(sum(abs(testimage - meanA)))];

end%for

for i=101:129,

filename = strcat('c',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhandswater,filename),'jpeg'));

P = [P; sum(sum(abs(testimage - meanA)))];

end%for

%%% Plot: Figure 2. L1 norm vs. the image index (ie. plot(1:length(P),P) ) and

%%% label all axis and title clearly

P = [Pa; Pb; P];

figure(2)

plot(P)

xlabel('Image Number')

ylabel('Pseudo Power')

title('Figure 2. Pseudo power for Faucet Controller Design')

grid on

%%% Q13. Above and below what power threshold values would you set to determine:

%%% a. when to turn on the water?

%%% b. when to turn off the water?

%%% Refer to Figure 2 and label these thresholds on the plot.

%%% Q14. Is the power a good indicator for our controller?

%%% What could we improve upon? Explain.

disp('part C')

% ************************************************************ %

% %

% Part C: Using image power and our masks to determine control %

% %

% ************************************************************ %

% The technique used in Part B could be improved upon using the masks

% we constructed n Part A. First of all, when we compute the image power

% we sum over all the pixels, and so the more noise in an image the more

% variable our power estimate. Second of all, we sometimes have difficulty

% determining when we are in a state of only water and only hands. It

% would be helpful if there was a way to know when hands were present but

% no water was flowing in the maintenance of the faucet.

% To get specific regional information we use the masks to isolate the

% useful location and then compute the power in those regions. So consider

% the design of a simple control, using the power level and masks to determine

% when to turn the water on and when to turn it off.

PAH = sum(sum(meanA .* maskHands));

PAW = sum(sum(meanA .* maskWater));

L=58;

PHa = NaN*ones(L,1); PW = NaN*ones(L,1);

PWa = NaN*ones(L,1); PW = NaN*ones(L,1);

for i=0:L-1,

filename = strcat('a',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirsink,filename),'jpeg'));

PHa(i+1) = sum(sum(abs(testimage - meanA).*maskHands))/PAH;

PWa(i+1) = sum(sum(abs(testimage - meanA).*maskWater))/PAW;

end%for

L=155;

PHb = NaN*ones(L,1); PW = NaN*ones(L,1);

PWb = NaN*ones(L,1); PW = NaN*ones(L,1);

for i=0:154,

filename = strcat('b',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirwater,filename),'jpeg'));

PHb(i+1) = sum(sum(abs(testimage - meanA).*maskHands))/PAH;

PWb(i+1) = sum(sum(abs(testimage - meanA).*maskWater))/PAW;

end%for

PH = [];

PW = [];

for i=24:41,

filename = strcat('d',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhands,filename),'jpeg'));

PH = [PH; sum(sum(abs(testimage - meanA).*maskHands))/PAH];

PW = [PW; sum(sum(abs(testimage - meanA).*maskWater))/PAW];

end%for

for i=101:129,

filename = strcat('c',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhandswater,filename),'jpeg'));

PH = [PH; sum(sum(abs(testimage - meanA).*maskHands))/PAH];

PW = [PW; sum(sum(abs(testimage - meanA).*maskWater))/PAW];

end%for

%%% 3x1 Plot: Figure 3. Label all axis and titles clearly

%%% 3a. [3,1] Power of image vs. image index

%%% 3b. [3,2] Normalized power of water-only region vs. image index

%%% 3c. [3,3] Normalized power of hands-only region vs. image index

PH = [PHa; PHb; PH];

PW = [PWa; PWb; PW];

figure(3)

subplot(3,1,1)

plot(P);grid on

ylabel('Pseudo Power');

title('Figure 3.1. Pseudo power for Faucet Controller Design')

subplot(3,1,2)

plot(PW);grid on

ylabel('Pseudo Power');

title('Figure 3.2. Normalized Pseudo power for Water Mask')

subplot(3,1,3)

plot(PH);grid on

xlabel('Image Number');

ylabel('Pseudo Power');

title('Figure 3.3. Normalized Pseudo power for Hand Mask')

%%% Q15. Why do we devide the power within the masks by the following:

%%% PAH = sum(sum(meanA .* maskHands));

%%% PAW = sum(sum(meanA .* maskWater));

%%% Q16. What logic would you use in your controller design now?

%%% Q17. Is there any benefit of using this approach to design over that in part B?

disp('part D')

% ******************************************************** %

% %

% Part D: Using probability and masks to determine control %

% %

% ******************************************************** %

% An even better approach than that given in Part B and C is to take

% advantage of probability theory and the standard deviations we computed

% earlier. The advantage to such an approach is that the control is based

% on a probabilistic interpretation of the image -- all indicators are

% between zero and one and the indicators have meaning. For instance, we could

% ask, what is the probability that the water is on, or that the hands are

% in the field of view.

% The code below is a simple attempt at design using a probability measure.

% PrH represents the probability that a hand is detected and PrW is the

% probability that water is detected.

PrH = [];

PrW = [];

countH = sum(sum(maskHands));

countW = sum(sum(maskWater));

for i=0:57,

filename = strcat('a',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirsink,filename),'jpeg'));

PrH = [PrH; sum(sum((abs(testimage - meanC) < 3*stdC) .* maskHands))/countH];

PrW = [PrW; sum(sum((abs(testimage - meanB) < 3*stdB) .* maskWater))/countW];

end%for

for i=0:154,

filename = strcat('b',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirwater,filename),'jpeg'));

PrH = [PrH; sum(sum((abs(testimage - meanC) < 3*stdC) .* maskHands))/countH];

PrW = [PrW; sum(sum((abs(testimage - meanB) < 3*stdB) .* maskWater))/countW];

end%for

for i=24:41,

filename = strcat('d',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhands,filename),'jpeg'));

PrH = [PrH; sum(sum((abs(testimage - meanC) < 3*stdC) .* maskHands))/countH];

PrW = [PrW; sum(sum((abs(testimage - meanB) < 3*stdB) .* maskWater))/countW];

end%for

for i=101:129,

filename = strcat('c',sprintf('%03d',i),'.',filetype);

testimage = double(imread(strcat(rootdir,filedirhandswater,filename),'jpeg'));

PrH = [PrH; sum(sum((abs(testimage - meanC) < 3*stdC) .* maskHands))/countH];

PrW = [PrW; sum(sum((abs(testimage - meanB) < 3*stdB) .* maskWater))/countW];

end%for

%%% 3x1 Plot: Figure 4. Label all axis and titles clearly

%%% 3a. [3,1] Power of image vs. image index (Part B graph)

%%% 3b. [3,2] Probability of detecting water [%] vs. image index

%%% 3c. [3,3] Probability of detecting hands [%] vs. image index

figure(4)

subplot(3,1,1)

plot(P);grid on

ylabel('Pseudo Power');

title('Figure 4.1. Pseudo power for Faucet Controller Design')

subplot(3,1,2)

plot(PrW);grid on

ylabel('Probability of Water');

title('Figure 4.2. Probability of Water in Image')

subplot(3,1,3)

plot(PrH);grid on

xlabel('Image Number');

ylabel('Probability of Hands');

title('Figure 4.3. Probability of Hands in Image')

%%% Q18. Explain the rational and logic behind the probability measures

%%% defined as:

%%% PrH = sum(sum((abs(testimage - meanC) < 3*stdC) .* maskHands))/countH;

%%% PrW = sum(sum((abs(testimage - meanB) < 3*stdB) .* maskWater))/countW;

%%% Q19. What logic would you use in your controller design now?

%%% Q20. Can you see any advatage in using this approach over that in part B

%%% and Part C?