% *************************************************************************
% ECE431F Digital Signal Processing
%
% Lab 4: Practical application of DSP to a real world problem:
% Processing digital signals from an infrared sensor operated faucet.
%
% Octave m. file scripts, you can rename it to lab4.m if you like.
%
% Please rename "imread_octave.m" to "imread.m"
%
% 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;
% 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";
% *********************************** %
% %
% 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 = loadjpg(strcat(rootdir,filedirsink,filename));
fprintf(stderr,"a mean %d\r",i);
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 = loadjpg(strcat(rootdir,filedirsink,filename));
fprintf(stderr,"a var %d\r",i);
varA = varA + (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);
disp('b mean calculated');
L=155; % number of measurements in the dataset
for i=0:L-1,
filename = strcat('b',sprintf('%03d',i),'.',filetype);
testimage = loadjpg(strcat(rootdir,filedirwater,filename));
fprintf(stderr,"b mean %d\r",i);
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),filetype));
testimage = loadjpg(strcat(rootdir,filedirwater,filename));
fprintf(stderr,"b var %d\r",i);
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 = loadjpg(strcat(rootdir,filedirhands,filename));
fprintf(stderr,"c mean %d\r",i);
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 = loadjpg(strcat(rootdir,filedirhands,filename));
fprintf(stderr,"c var %d\r",i);
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
imagesc(meanA,1);colormap(gray);title('1.1. The Sink')
imagesc(stdA,1);colormap(gray);title('1.2. Standard Deviation of Sink')
imagesc(maskHands,1);colormap(gray);title('1.3. Mask of Hands Only')
imagesc(meanB,1);colormap(gray);title('1.4. The Water')
imagesc(stdB,1);colormap(gray);title('1.5. Standard Deviation of Water')
imagesc(maskWater,1);colormap(gray);title('1.6. Mask of Water Only')
imagesc(meanC,1);colormap(gray);title('1.7. The Hands')
imagesc(stdC,1);colormap(gray);title('1.8. Standard Deviation of Hands')
imagesc(maskC,1);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 = loadjpg(strcat(rootdir,filedirsink,filename));
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 = loadjpg(strcat(rootdir,filedirwater,filename));
Pb(i+1) = sum(sum(abs(testimage - meanA)));
end%for
P =[];
for i=24:41,
filename = strcat('d',sprintf('%03d',i),'.',filetype);
testimage = loadjpg(strcat(rootdir,filedirhands,filename));
P = [P; sum(sum(abs(testimage - meanA)))];
end%for
for i=101:129,
filename = strcat('c',sprintf('%03d',i),'.',filetype);
testimage = loadjpg(strcat(rootdir,filedirhandswater,filename));
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)
xlabel('Image Number')
ylabel('Pseudo Power')
title('Figure 2. Pseudo power for Faucet Controller Design')
plot(P)
%%% 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 = loadjpg(strcat(rootdir,filedirsink,filename));
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 = loadjpg(strcat(rootdir,filedirwater,filename));
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 = loadjpg(strcat(rootdir,filedirhands,filename));
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 = loadjpg(strcat(rootdir,filedirhandswater,filename));
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)
ylabel('Pseudo Power');
title('Figure 3.1. Pseudo power for Faucet Controller Design')
plot(P);
subplot(3,1,2)
ylabel('Pseudo Power');
title('Figure 3.2. Normalized Pseudo power for Water Mask')
plot(PW);
subplot(3,1,3)
xlabel('Image Number');
ylabel('Pseudo Power');
title('Figure 3.3. Normalized Pseudo power for Hand Mask')
plot(PH);
%%% 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 = loadjpg(strcat(rootdir,filedirsink,filename));
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 = loadjpg(strcat(rootdir,filedirwater,filename));
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 = loadjpg(strcat(rootdir,filedirhands,filename));
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 = loadjpg(strcat(rootdir,filedirhandswater,filename));
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)
ylabel('Pseudo Power');
title('Figure 4.1. Pseudo power for Faucet Controller Design')
plot(P);
subplot(3,1,2)
ylabel('Probability of Water');
title('Figure 4.2. Probability of Water in Image')
plot(PrW);
subplot(3,1,3)
xlabel('Image Number');
ylabel('Probability of Hands');
title('Figure 4.3. Probability of Hands in Image')
plot(PrH);
%%% 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?