Commit after merge

Author: iuya

Motion-Estimation/Block_Matching.m

@@ -0,0 +1,51 @@
+function [ x_predicted, y_predicted, min_MSE ] = Block_Matching( reference_frame, new_frame, p, x, y, block_size  )
+% Block_Matching Algorithm. Block Matching between 2 subsequent frames to detect motion estimation
+%
+% Input
+%   reference_frame :  The frame{t-1}
+%   new_frame       :  The frame{t}
+%   p               :  Search parameter
+%   x,y             :  Position withing the frame
+%   block_size      :  Size of the macroblock
+% 
+% Ouput 
+%   motion_vect     :  Motion vector for each macroblock
+%   min_MSE         :  The minimum MSE (Mean Square Error)
+% 
+
+reference_block = reference_frame(x:x+(block_size-1),y:y+(block_size-1));
+MSE = ones(2*p + 1, 2*p +1) * 65537;
+
+% Search the correspondence of the macroblock of a reference frame
+% using an exhaustive method. This method searches for the minimum
+% cost function at each possible location with respect to the search parameter.
+
+for m=-p:p
+    for n=-p:p
+        new_xpos = x + m;
+        new_ypos = y + n;
+        if(new_xpos+block_size>size(reference_frame,1) || new_ypos+block_size>size(reference_frame,2) || new_xpos<1 || new_ypos<1)
+            continue;
+        else
+            new_block =new_frame(new_xpos:new_xpos+(block_size-1),new_ypos:new_ypos+(block_size-1));
+            MSE(m+p+1,n+p+1)=double(mean2(((reference_block-new_block).^2)));
+        end
+    end
+end
+
+% Create the motion vectors with the minimum cost function.
+min = 70000;
+for i = 1:size(MSE,1)
+    for j = 1:size(MSE,2)
+        if (MSE(i,j) < min)
+            min = MSE(i,j);dx = j; dy = i;                    
+        end
+    end
+end
+
+min_MSE = min;
+x_predicted = dx-p-1;
+y_predicted = dy-p-1;
+
+end
+

Motion-Estimation/Motion_Compensation.m

@@ -0,0 +1,34 @@
+function  estimated_image = Motion_Compensation( reference_frame, motion_vect, block_size)
+%Motion_Compensation algorithm. Predicts the next image using the motion
+%vectors
+%   
+%  Input 
+%   reference_frame  : The frame{t-1} 
+%   motion_vect      : The motion vectors for each macroblock 
+%   block_size       : Size of the macroblock
+%
+%  Output
+%   estimated_image  : Result of using the motion vectors to predict a frame in a video stream
+%
+
+[row, col] = size(reference_frame);
+vect_count = 1;
+image_compensation = reference_frame;
+
+for i = 1:block_size:row-block_size
+    for j = 1:block_size:col-block_size        
+        dy = motion_vect(1,vect_count);
+        dx = motion_vect(2,vect_count);
+        y_block = i + dy;
+        x_block = j + dx;
+        if(y_block>=1&&x_block>=1&&dy~=0&&dx~=0)      
+            image_compensation(i:i+block_size-1,j:j+block_size-1) = reference_frame(y_block:y_block+block_size-1, x_block:x_block+block_size-1);
+        end
+        vect_count = vect_count + 1;
+    end
+end
+
+estimated_image = image_compensation;
+
+end
+

Motion-Estimation/Motion_Estimation.m

@@ -0,0 +1,28 @@
+function motion_vectors = Motion_Estimation(reference_frame, new_frame, block_size, p)
+%Motion_Estimation algorithm  Computes motion vectors using exhaustive search method
+%
+% Input
+%   reference_frame : The frame{t-1} 
+%   new_frame       : The frame{t}
+%   p               : Search parameter
+%   block_size      : Size of the macroblock
+%
+% Ouput
+%   motion_vectors  : The motion vectors for each macroblock 
+%
+
+[rows, cols] = size(reference_frame);
+motion_vectors = zeros(3,rows*cols/block_size^2);
+
+vect_count = 1;
+for posx = 1:block_size:rows-block_size
+    for posy = 1:block_size:cols-block_size
+          [predicted_posx,predicted_posy,min_MSE]=Block_Matching( reference_frame, new_frame, p, posx, posy, block_size  );
+          motion_vectors(1:2,vect_count) = [predicted_posx,predicted_posy];
+          motion_vectors(3,vect_count) = min_MSE;
+          vect_count = vect_count+1;
+    end
+end
+
+end
+

Motion-Estimation/Tests/40-41.svg

@@ -0,0 +1,23 @@
+function draw_Optical_Flow( reference_frame, motion_vect, block_size )
+%UNTITLED Draw the motion vectors in the reference frame
+
+[row, col] = size(reference_frame);
+mbCount = 1;
+subplot(2,2,4),subimage(reference_frame),title('Motion Vectors');
+hold on;
+for i = 1:block_size:row-block_size
+    for j = 1:block_size:col-block_size       
+        dy = motion_vect(1,mbCount);
+        dx = motion_vect(2,mbCount);
+        y_block = i + dy;
+        x_block = j + dx;
+        if((y_block>=1&&x_block>=1&&dy~=0&&dx~=0) &&  (motion_vect(3,mbCount)>0.4) && (x_block<=size(reference_frame,1)) && (y_block<=size(reference_frame,2)))
+            quiver(j,i,dx, dy,'color', 'b', 'linewidth', 1);
+        end
+        mbCount = mbCount + 1;
+    end
+end
+hold off;
+
+end
+

Motion-Estimation/draw_Optical_Flow.m

@@ -0,0 +1,11 @@
+function PSNR = image_PSNR(reference_image, estimated_image)
+
+[row, col] = size(reference_image);
+sum_err = 0;
+for i = 1:row
+    for j = 1:col
+        sum_err = sum_err + (reference_image(i,j) - estimated_image(i,j))^2;
+    end
+end
+mse = sum_err/(row*col);
+PSNR = 10*log10(double(255^2/mse));

Motion-Estimation/image_PSNR.m

@@ -0,0 +1,39 @@
+% Read the video file
+video=VideoReader('../Test-images/Xylophone.mp4');
+n=video.NumberofFrames;
+
+% Divide the video in frames
+for x=1:n
+frame=read(video,x);
+imwrite(frame,sprintf('../Test-images/Frames/image%d.jpg',x));
+end
+
+% Compute the motion vectors between 2 frames
+reference_frame = imread('../Test-images/Frames/image40.jpg');
+new_frame = imread('../Test-images/Frames/image41.jpg');
+
+% Transform in gray scale for simplification
+reference_frame = rgb2gray(reference_frame);
+new_frame = rgb2gray(new_frame);
+
+% Plot the reference and new frames
+subplot(2,2,1),subimage(reference_frame),title('Reference Frame'); hold on;
+subplot(2,2,2),subimage(new_frame),title('New Frame');
+
+% p is the search parameter, block_size is the size of the macroblock
+p=7;block_size=16;
+
+% Motion estimation using an exhaustive method
+[motion_vect] = Motion_Estimation(reference_frame,new_frame,block_size,p);
+
+% Motion competation using the motion vectors 
+estimated_frame= Motion_Compensation(reference_frame, motion_vect, block_size);
+subplot(2,2,3),subimage(estimated_frame),title('Estimated Frame');
+draw_Optical_Flow( reference_frame, motion_vect, block_size );
+hold off;
+
+% Peak Signal-to-Noise Ratio (PSNR) to measure the distortion of the frame which is estimated and the predicted one
+PSNR = image_PSNR(new_frame, estimated_frame);
+fprintf('Peak-Signal-To-Noise Ratio: %d \n',PSNR);
+
+

Motion-Estimation/script.m

@@ -0,0 +1,7 @@
+function paddedMatrix = paddingZeros(matrix, divisor)
+    %paddedMatrix = zeros(size(matrix,1), size(matrix,2), size(matrix,3));
+    rows = size(matrix,1);
+    cols = size(matrix,2);
+    paddedMatrix = wextend('ar','sym',matrix, round(rows/divisor)*divisor - rows, 'd');
+    paddedMatrix = wextend('ac','sym',paddedMatrix, round(cols/divisor)*divisor - cols, 'l');
+end

Stereoscopy/paddingZeros.m

@@ -0,0 +1,183 @@
+[img1, colormap] = imread('../Test-images/stereoscopy-images/image3.jpg');
+figure; image(img1), title('Original Image'); % Original image
+
+%  Split them into 2 separate images
+leftImage = img1(:,1:(end/2),:);
+
+rightImage = img1(:,end/2+1:end,:);
+figure; image(leftImage), title('Left Image');
+figure; image(rightImage), title('Right Image');
+
+leftImage = paddingZeros(leftImage,4); %We need a dimension divisible by 4
+rightImage = paddingZeros(rightImage,4);
+
+
+% Encode each eye's image using filters of different (usually chromatically opposite) colors, red and blue
+r = zeros(size(leftImage));
+gb = zeros(size(rightImage));
+r(:,:,1) = double(leftImage(:,:,1));
+gb(:,:,2:3) = double(rightImage(:,:,2:3));
+anaglyph = uint8(r+gb);
+figure,image(uint8(r)), title('Red')
+figure,image(uint8(gb)), title('Cyan')
+figure,image(anaglyph), title('My Anaglyph')
+
+% Create stereo anaglyph with Matlab function to check the result
+J = stereoAnaglyph(leftImage,rightImage);
+figure, image(uint8(J)), title('Anaglyph using Matlab function')
+
+% Use TD3 to code them using 2 decomposition layers
+leftImage_haar = zeros(size(r,1),size(r,2),size(r,3));
+rightImage_haar =zeros(size(gb,1),size(gb,2),size(gb,3));
+rows_l=size(leftImage,1);
+rows_r=size(rightImage,1);
+cols_l=size(leftImage,2);
+cols_r=size(rightImage,2);
+sub11_l = zeros(rows_l/4,cols_l/4,3); sub11_r = zeros(rows_r/4,cols_r/4,3);
+sub12_l = zeros(rows_l/4,cols_l/4,3); sub12_r = zeros(rows_r/4,cols_r/4,3);
+sub13_l = zeros(rows_l/4,cols_l/4,3); sub13_r = zeros(rows_r/4,cols_r/4,3);
+sub14_l = zeros(rows_l/4,cols_l/4,3); sub14_r = zeros(rows_r/4,cols_r/4,3);
+sub21_l = zeros(rows_l/2,cols_l/2,3); sub21_r = zeros(rows_r/2,cols_r/2,3);
+sub22_l = zeros(rows_l/2,cols_l/2,3); sub22_r = zeros(rows_r/2,cols_r/2,3);
+sub23_l = zeros(rows_l/2,cols_l/2,3); sub23_r = zeros(rows_r/2,cols_r/2,3);
+
+for i=1:3
+   % temp = leftImage(:,:,i);
+    temp = r(:,:,i);
+
+    leftImage_haar(:,:,i) = haar_transform_multilevel(temp,2);
+
+    sub11_l(:,:,i) = leftImage_haar(1:rows_l/4, 1:cols_l/4,i);
+    sub12_l(:,:,i) = leftImage_haar(1:rows_l/4, cols_l/4+1:cols_l/2,i);
+    sub13_l(:,:,i) = leftImage_haar(rows_l/4+1:rows_l/2, 1:cols_l/4,i);
+    sub14_l(:,:,i) = leftImage_haar(rows_l/4+1:rows_l/2, cols_l/4+1:cols_l/2,i);
+    sub21_l(:,:,i) = leftImage_haar(1:rows_l/2, cols_l/2+1:cols_l,i);
+    sub22_l(:,:,i) = leftImage_haar(rows_l/2+1:rows_l, 1:cols_l/2,i);
+    sub23_l(:,:,i) = leftImage_haar(rows_l/2+1:rows_l, cols_l/2+1:cols_l,i);
+
+    temp2 = gb(:,:,i);
+    rightImage_haar(:,:,i) = haar_transform_multilevel(temp2,2);
+    sub11_r(:,:,i) = rightImage_haar(1:rows_l/4, 1:cols_l/4,i);
+    sub12_r(:,:,i) = rightImage_haar(1:rows_l/4, cols_l/4+1:cols_l/2,i);
+    sub13_r(:,:,i) = rightImage_haar(rows_l/4+1:rows_l/2, 1:cols_l/4,i);
+    sub14_r(:,:,i) = rightImage_haar(rows_l/4+1:rows_l/2, cols_l/4+1:cols_l/2,i);
+    sub21_r(:,:,i) = rightImage_haar(1:rows_l/2, cols_l/2+1:cols_l,i);
+    sub22_r(:,:,i) = rightImage_haar(rows_l/2+1:rows_l, 1:cols_l/2,i);
+    sub23_r(:,:,i) = rightImage_haar(rows_l/2+1:rows_l, cols_l/2+1:cols_l,i);
+end
+figure;image(uint8(leftImage_haar));title('Left Image after haar filter');
+figure;image(uint8(rightImage_haar));title('Right Image after haar filter');
+
+% With more detail
+figure
+subplot(4,4,1)
+image(sub11_l./255)
+subplot(4,4,2)
+image(sub12_l)
+subplot(4,4,5)
+image(sub13_l)
+subplot(4,4,6)
+image(sub14_l)
+% Now for the bigger subbands we are going to use bigger subplots
+subplot(4,4,[3:4 7:8])
+image(sub21_l)
+subplot(4,4,[9:10 13:14])
+image(sub22_l)
+subplot(4,4,[11:12 15:16])
+image(sub23_l)
+
+% With more detail
+figure
+subplot(4,4,1)
+image(sub11_r./255)
+subplot(4,4,2)
+image(sub12_r)
+subplot(4,4,5)
+image(sub13_r)
+subplot(4,4,6)
+image(sub14_r)
+% Now for the bigger subbands we are going to use bigger subplots
+subplot(4,4,[3:4 7:8])
+image(sub21_r)
+subplot(4,4,[9:10 13:14])
+image(sub22_r)
+subplot(4,4,[11:12 15:16])
+image(sub23_r)
+
+leftImage_dq = zeros(size(leftImage,1),size(leftImage,2),size(leftImage,3));
+leftImage_synth = zeros(size(leftImage,1),size(leftImage,2),size(leftImage,3));
+
+rightImage_dq =zeros(size(rightImage,1),size(rightImage,2),size(rightImage,3));
+rightImage_synth = zeros(size(rightImage,1),size(rightImage,2),size(rightImage,3));
+
+for i=1:3
+    q_sub11_l(:,:,i) = sub11_l(:,:,i);
+    q_sub12_l(:,:,i) = quantize_matrix(sub12_l(:,:,i),4);
+    q_sub13_l(:,:,i) = quantize_matrix(sub13_l(:,:,i),4);
+    q_sub14_l(:,:,i) = quantize_matrix(sub14_l(:,:,i),4);
+    q_sub21_l(:,:,i) = quantize_matrix(sub21_l(:,:,i),3);
+    q_sub22_l(:,:,i) = quantize_matrix(sub22_l(:,:,i),3);
+    q_sub23_l(:,:,i) = quantize_matrix(sub23_l(:,:,i),3);
+
+    q_sub11_r(:,:,i) = sub11_r(:,:,i);
+    q_sub12_r(:,:,i) = quantize_matrix(sub12_r(:,:,i),4);
+    q_sub13_r(:,:,i) = quantize_matrix(sub13_r(:,:,i),4);
+    q_sub14_r(:,:,i) = quantize_matrix(sub14_r(:,:,i),4);
+    q_sub21_r(:,:,i) = quantize_matrix(sub21_r(:,:,i),3);
+    q_sub22_r(:,:,i) = quantize_matrix(sub22_r(:,:,i),3);
+    q_sub23_r(:,:,i) = quantize_matrix(sub23_r(:,:,i),3);
+    
+    %Synthesis
+    leftImage_dq(1:rows_l/4, 1:cols_l/4,i)= q_sub11_l(:,:,i);
+    leftImage_dq(1:rows_l/4, cols_l/4+1:cols_l/2,i) = dequantize_matrix(q_sub12_l(:,:,i),4,8);
+    leftImage_dq(rows_l/4+1:rows_l/2, 1:cols_l/4,i) = dequantize_matrix(q_sub13_l(:,:,i),4,8);
+    leftImage_dq(rows_l/4+1:rows_l/2, cols_l/4+1:cols_l/2,i) = dequantize_matrix(q_sub14_l(:,:,i),4,8);
+    leftImage_dq(1:rows_l/2, cols_l/2+1:cols_l,i) = dequantize_matrix(q_sub21_l(:,:,i),3,8);
+    leftImage_dq(rows_l/2+1:rows_l, 1:cols_l/2,i) = dequantize_matrix(q_sub22_l(:,:,i),3,8);
+    leftImage_dq(rows_l/2+1:rows_l, cols_l/2+1:cols_l,i) = dequantize_matrix(q_sub23_l(:,:,i),4,8);
+    leftImage_synth(:,:,i) = haar_reverse_multilevel(leftImage_dq(:,:,i),2);
+
+    rightImage_dq(1:rows_r/4, 1:cols_r/4,i)= q_sub11_r(:,:,i);
+    rightImage_dq(1:rows_r/4, cols_r/4+1:cols_r/2,i) = dequantize_matrix(q_sub12_r(:,:,i),4,8);
+    rightImage_dq(rows_r/4+1:rows_r/2, 1:cols_r/4,i) = dequantize_matrix(q_sub13_r(:,:,i),4,8);
+    rightImage_dq(rows_r/4+1:rows_r/2, cols_r/4+1:cols_r/2,i) = dequantize_matrix(q_sub14_r(:,:,i),4,8);
+    rightImage_dq(1:rows_r/2, cols_r/2+1:cols_r,i) = dequantize_matrix(q_sub21_r(:,:,i),3,8);
+    rightImage_dq(rows_r/2+1:rows_r, 1:cols_r/2,i) = dequantize_matrix(q_sub22_r(:,:,i),3,8);
+    rightImage_dq(rows_r/2+1:rows_r, cols_r/2+1:cols_r,i) = dequantize_matrix(q_sub23_r(:,:,i),4,8);
+    rightImage_synth(:,:,i) = haar_reverse_multilevel(rightImage_dq(:,:,i),2);
+end
+<<<<<<< HEAD
+figure; image(uint8(rightImage_synth));
+figure; image(uint8(leftImage_synth));
+
+%Anaglyph of the synthesis
+anaglyph_synth = uint8(leftImage_synth+rightImage_synth);
+figure,image(anaglyph_synth), title('Synthesized Anaglyph')
+
+function paddedMatrix = paddingZeros(matrix, divisor)
+    %paddedMatrix = zeros(size(matrix,1), size(matrix,2), size(matrix,3));
+    rows = size(matrix,1);
+    cols = size(matrix,2);
+    paddedMatrix = wextend('ar','sym',matrix, round(rows/divisor)*divisor - rows, 'd');
+    paddedMatrix = wextend('ac','sym',paddedMatrix, round(cols/divisor)*divisor - cols, 'l');
+end
+=======
+figure; image(uint8(rightImage_synth)), title('Right Image Synth');
+figure; image(uint8(leftImage_synth)), title('Left Image Synth');
+
+% Encode each eye's image using filters of different (usually chromatically opposite) colors, red and blue
+rightImage_synth  = imtranslate(rightImage_synth,[10, 0]);
+
+figure,image(rightImage_synth), title('Anaglyph Synth');
+r = zeros(size(rightImage_synth));
+gb = zeros(size(leftImage_synth));
+r(:,:,1) = double(rightImage_synth(:,:,1));
+gb(:,:,2:3) = double(leftImage_synth(:,:,2:3));
+anaglyph = uint8(r+gb);
+figure,image(anaglyph), title('Anaglyph Synth');
+
+% Create stereo anaglyph with Matlab function to check the result
+J = stereoAnaglyph(rightImage_synth,leftImage_synth);
+figure, image(uint8(J)), title('Anaglyph using Matlab function')
+
+>>>>>>> 6a9213beae2e0128673f00ecc9ab774415114025