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

% 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?