# Laborator 3 SCCV

https://base-n.de/matlab/code_beautifier.html

https://typora.io/
https://www.sublimetext.com/

https://www.markdownguide.org/getting-started/
https://www.markdownguide.org/basic-syntax/

### Aplicatia 1 Lab3 SCCV

> Radu - cod netestat
(Melania - testat)
Iris - modificat acolo la C(m+1,n+1) cu cos in loc de sin

clear, clc, close all;
for N = [4 8]
    C = zeros(N, N); % N este dimensiunea blocului transformatei DCT
    % Pentru cazul prezent blocurile sunte de dimensiune 4x4 si 8x8
    % Se calculează coeficienţii C
    for m = 0:1:N-1
        for n = 0:1:N-1
            if n == 0
                k = sqrt (1 / N);
            else
                k = sqrt (2 / N);
            end
            C(m+1, n+1) = k * cos(((2*m+1)*n*pi) / (N+1));
        end
    end
    % Se determina funcţiile de baza
    figure;
    colormap('gray');
    for m = 0:1:N-1
        for n = 0:1:N-1
            subplot(N, N, m * N + n + 1);
            Y = [zeros(m, N);
            zeros(1,n) 1 zeros(1, N - n - 1);
            zeros(N - m - 1, N)];
            X = C * Y * C';
            imagesc (X);
            axis square;
            axis off;
        end
    end
end

### Exercitiul 1 Lab 3 SCCV
Transformata sin
(Iris)

clear, clc, close all;
for N = [4 8]
    C = zeros(N, N); % N este dimensiunea blocului transformatei DCT
    % Pentru cazul prezent blocurile sunte de dimensiune 4x4 si 8x8
    % Se calculează coeficienţii C
    for m = 0:1:N-1
        for n = 0:1:N-1
                k = sqrt (2 / N+1);
            C(m+1, n+1) = k * sin(((m+1)*(n+1)*pi) / (N+1));
        end
    end
    % Se determina funcţiile de baza
    figure;
    colormap('gray');
    for m = 0:1:N-1
        for n = 0:1:N-1
            subplot(N, N, m * N + n + 1);
            Y = [zeros(m, N);
            zeros(1,n) 1 zeros(1, N - n - 1);
            zeros(N - m - 1, N)];
            X = C * Y * C';
            imagesc (X);
            axis square;
            axis off;
        end
    end
end

### Aplicatia 2 Lab3 SCCV

Pentru 2D-FFT:
(Melania - testat)

clear, clc, close all;
x = imread('baboon200_200.jpg');
y = rgb2gray(x);
[inaltime, latime] = size(y);
Coef = fft2 (y);
imshow(log(abs(fftshift(Coef)) + 1), []), colormap(jet (16)), colorbar

Pentru 2D-DCT:
(Melania - testat)

x = imread('baboon200_200.jpg');
y = rgb2gray(x);
[inaltime, latime] = size(y);
Coef = dct2(y);
imshow(log(abs(Coef)), []), colormap(jet(16)), colorbar

### Aplicatia 3 Lab3 SCCV

> Radu
> (Melania - testat)

function G = TCgain(x, N)
    % G = TCgain(x)
    % Se calculeaza castigul codarii cand se aplica 2D-DCT
    % Variabile de intrare:
    % x = Imaginea de intensitate
    % N = dimensiunea blocului de transformare (de regula patratic)
    % Se obtine la iesire:
    % G = câştigul de codare in dB raportat la PCM
    [Height, Width] = size(x);
    fun = @dct2;
    Y = blkproc(x, [N N], fun); % Se calc 2D-DCT de dimensiune N x N
    % Se calculează varianta coeficienţilor DCT
    CoefVar = zeros(N, N);
    for i = 1:N
        for j = 1:N
            t = std2(Y(i:N:Height, j:N:Width));
            CoefVar(i, j) = t * t;
        end
    end
    
    MeanVar = mean(CoefVar(:));%Media variantei
    P1 = CoefVar.^(1/(N^2));
    GeoMeanVar = prod(P1(:)); % Media geometrica a variantei coeficienţilor
    G = 10*log10(MeanVar/GeoMeanVar); % Câştigul
end

Apelarea functiei:
(Melania - testat)

clear, clc, close all;
A = imread('baboon.jpg');
% imshow(x)
[Height, Width, Depth] = size(A);
if Depth == 1
    A1 = double(A);
else
    A = double(rgb2ycbcr(A));
    A1 = A(:, :, 1);
end

% figure, imshow(A1, []), title('Baboon');

N = 4;
for j = 1:4
    G = TCgain(A1, N);
    sprintf('Castigul codarii pentru %dx%d DCT = %4.2f', N, N, G)
    N = 2*N
end

### Aplicatia 4 Lab 3 SCCV
> Denisa - cod netestat
> Radu - testat, rescris/corectat

clear, clc, close all;
A = imread('..\PozePrelucrare\lena-256x256.jpg');
imshow(A)
[Height,Width,Depth] = size(A);
N = 8; % dimensiunea blocului
if mod(Height, N) ~= 0
    Height = floor(Height/N) * N;
end
if mod(Width, N) ~= 0
    Width = floor(Width/N) * N ;
end

Al = A(1:Height, 1:Width, :);

clear A
A = Al;

SamplingFormat = '4:2:0';
if Depth == 1
    y = double(A);
else
    %Conversie din RGB in YCbCr
    A = double(rgb2ycbcr(A));
    y = A(:,:,1);

    % Sub-esantionarea imformatiei de crominanta
    switch SamplingFormat
        case '4:2:0'
            Cb = imresize(A(:,:,2),[Height/2 Width/2], 'bicubic');
            Cr = imresize(A(:,:,3),[Height/2 Width/2], 'bicubic');
        case '4:2:2' 
            Cb = imresize(A(:,1i,2),[Height Width/2], 'bicubic');
            Cr = imresize(A(:,:,3),[Height Width/2], 'bicubic');
    end
end


% Matrice de cuantizare informaţie luminanta
jpgQstepsY = [16 11 10 16 24 40 51 61;...
              12 12 14 19 26 58 60 55;...
              14 13 16 24 40 57 69 56;...
              14 17 22 29 51 87 80 62;...
              18 22 37 56 68 109 103 77;...
              24 35 55 64 81 104 113 92;...
              49 64 78 87 103 121 120 101;...
              72 92 95 98 112 100 103 99];
          
QstepsY = jpgQstepsY;
Qscale = 1.5;
Yy = zeros(N,N);
xqY = zeros(Height,Width);
acBitsY = 0;
dcBitsY = 0;
          
if Depth > 1
    jpgQstepsC = [17 18 24 47 66 99 99 99;...
                  18 21 26 66 99 99 99 99;...
                  24 26 56 99 99 99 99 99;...
                  47 66 99 99 99 99 99 99;...
                  99 99 99 99 99 99 99 99;...
                  99 99 99 99 99 99 99 99;...
                  99 99 99 99 99 99 99 99;...
                  99 99 99 99 99 99 99 99];
    QstepsC = jpgQstepsC;
    YCb = zeros(N,N);
    YCr = zeros(N,N);
    switch SamplingFormat
        case '4:2:0'
            xqCb = zeros(Height/2, Width/2);
            xqCr = zeros(Height/2, Width/2);
        case '4:2:2'
            xqCb = zeros(Height, Width/2);
            xqCr = zeros(Height, Width/2);
    end
    acBitsCb = 0;
    dcBitsCb = 0;
    acBitsCr = 0;
    dcBitsCr = 0;
end

% 5e calculează cuvinetele de cod asociate fiecărei componente
% Pentru componenta de luminanta
for m = 1:N:Height
  for n = 1:N:Width
      t = y(m:m+N-1,n:n+N-1)-128;
      Yy = dct2(t); % Se aplica transformata DCT/2D pentru blocuri
                    % de dimensiune N x N
      % Se cuantizeaza coeficienţii DCT
      temp = floor(Yy ./ (Qscale*QstepsY) + 0.5);
      % Se calculează biţii de cod pentru diferentele DC
      if n == 1
         DC = temp(1,1);
         dcBitsY = dcBitsY + jpgDCbits(DC, 'Y');
      else
         DC = temp(1,1) - DC;
         dcBitsY = dcBitsY + jpgDCbits(DC, 'Y');
         DC = temp(1,1);
      end
      % Se calculează biţii de cod pentru coeficienţii AC
      ACblkBits = jpgACbits(temp,'Y');
      acBitsY = acBitsY + ACblkBits;
      % Se decuantizeaza coeficienţii si se calculează IDCT
      xqY(m:m+N-1,n:n+N-1)= idct2(temp.*(Qscale*QstepsY))+ 128;
  end
end

% Daca imaginea aplicata la intrare este color
% se calculează biţii de cod pentru componentele de crominanta
if Depth > 1
    if strcmpi(SamplingFormat,'4:2:0')
        EndRow = Height/2;
    else
        EndRow = Height;
    end
    for m = 1:N:EndRow
       for n = 1:N:Width/2
           tl = Cb(m:m+N-1,n:n+N-1) - 128;
           t2 = Cr(m:m+N-1,n:n+N-1) - 128;
           Ycb = dct2(tl); % Transformata DCT de dimensiune N x N
           Ycr = dct2(t2);
           temp1 = floor(Ycb ./ (Qscale*QstepsC) + 0.5);
           temp2 = floor(Ycr ./ (Qscale*QstepsC) + 0.5);
           if n == 1
              DC1 = temp1(1, 1);
              DC2 = temp2(1, 1);
              dcBitsCb = dcBitsCb + jpgDCbits(DC1, 'C');
              dcBitsCr = dcBitsCr + jpgDCbits(DC2, 'C');
           else
              DC1 = temp1(1,1) - DC1;
              DC2 = temp2(1,1) - DC2;
              dcBitsCb = dcBitsCb + jpgDCbits(DC1, 'C');
              dcBitsCr = dcBitsCr + jpgDCbits(DC2, 'C');
              DC1 = temp1(1,1);
              DC2 = temp2(1,1);
           end
           ACblkBitsl = jpgACbits(temp1, 'C');
           ACblkBits2 = jpgACbits(temp2, 'C');
           acBitsCb = acBitsCb + ACblkBitsl;
           acBitsCr = acBitsCr + ACblkBits2;
           
           % Decuantizarea coeficienţilor si aplicarea IDCT
           xqCb(m:m+N-1,n:n+N-1)= idct2(temp1 .* (Qscale*QstepsC))+ 128;
           xqCr(m:m+N-1,n:n+N-1)= idct2(temp2 .* (Qscale*QstepsC))+ 128;
       end
    end
end

mse = std2(y - xqY);
snr = 20 * log10(std2(y) / mse);
fprintf('SNR = %4.2f\n', snr)
if Depth == 1
    TotalBits = acBitsY + dcBitsY;
    figure, imshow(xqY, [])
    title(['JPG compressed ' '@ ' num2str(TotalBits/(Height*Width)) 'bpp'])
else
   TotalBits = acBitsY + dcBitsY + dcBitsCb +...
               acBitsCb + dcBitsCr + acBitsCr;
   cl = imresize(xqCb,[Height Width], 'bicubic');
   c2 = imresize(xqCr,[Height Width], 'bicubic');
   mseb = std2(A(:,:,2) - cl);
   snrb = 20*log10(std2(A(:,:,2)) / mseb);
   msec = std2(A(:,:,3) - c2);
   snrc = 20*log10(std2(A(:,:,3)) / msec);
   fprintf('SNR(Cb)=%4.2fdB\tSNR(Cr)=%4.2fdB\n',snrb,snrc)
   xq(:,:,1) = xqY;
   xq(:,:,2) = cl;
   xq(:,:,3) = c2;
   figure, imshow(ycbcr2rgb(uint8(round(xq))))
   title(['JPG compressed ' '@' num2str(TotalBits/(Height *Width)) ' bpp'])
end
fprintf('Bit rate = %4.2f bpp\n', TotalBits/(Height*Width))

Functii aplicatia 4
> Madalin - cod testat

function Bits = jpgDCbits(dc,C)
%Determina numarul exact de biti necesar codarii entropice pentru
%coeficientii DC ce au fost codati deferential in cursul codarii JPEG
%ce utilizeaza tabelele predefinite Huffman
%Variabile de intrare:
%dc - valorile prezise pt coef DC
%C = 'Y' pentru luminanta si 'C' pt crominanta
%Variabile la iesire:
%Nr de biti aferent codarii diferentiale
%CodeLengthY contine lung cuvintelor de cod pt luminanta DC diff pt categ
%de la 0 la 11 inclusiv
%CodeLengthC contine lung cuvintelor de cod pt crominanta DC

CodeLengthY = int16([3 4 5 5 7 8 10 12 14 16 18 20]);
CodeLengthC = int16([3 4 5 6 7 9 11 13 15 18 19 20]);
switch C
    case 'Y'
        CodeLength = CodeLengthY;
    case 'C'
        CodeLength = CodeLengthC;
end
if dc==0
    Bits = double(CodeLength(1));
else
    Bits = double(CodeLength(round(log2(abs(dc))+0.5)+1));
end

function Bits = jpgACbits(x,C)
RLCy = cell(1,15);
RLCy{1} = int16([4 3 4 6 8 10 12 14 18 25 26]);
RLCy{2} = int16([5 8 10 13 16 22 23 24 25 26]);
RLCy{3} = int16([6 10 13 20 21 22 23 24 25 26]);
RLCy{4} = int16([7 11 14 20 21 22 23 24 25 26]);
RLCy{5} = int16([7 12 19 20 21 22 23 24 25 26]);
RLCy{6} = int16([8 12 19 20 21 22 23 24 25 26]);
RLCy{7} = int16([8 13 19 20 21 22 23 24 25 26]);
RLCy{8} = int16([9 13 19 20 21 22 23 24 25 26]);
RLCy{9} = int16([9 17 19 20 21 22 23 24 25 26]);
RLCy{10} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{11} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{12} = int16([10 18 19 20 21 22 23 24 25 26]);
RLCy{13} = int16([11 18 19 20 21 22 23 24 25 26]);
RLCy{14} = int16([12 18 19 20 21 22 23 24 25 26]);
RLCy{15} = int16([13 18 19 20 21 22 23 24 25 26]);
RLCy{16} = int16([12 17 18 19 20 21 22 23 24 25 26]);

RLCc = cell(1,15);
RLCc{1} = int16([2 3 5 7 9 12 15 23 24 25 26]);
RLCc{2} = int16([5 8 11 14 20 22 23 24 25 26]);
RLCc{3} = int16([6 9 13 15 21 22 23 24 25 26]);
RLCc{4} = int16([6 9 12 15 21 22 23 24 25 26]);
RLCc{5} = int16([6 10 14 20 21 22 23 24 25 26]);
RLCc{6} = int16([7 11 17 20 21 22 23 24 25 26]);
RLCc{7} = int16([8 12 19 20 21 22 23 24 25 26]);
RLCc{8} = int16([7 12 19 20 21 22 23 24 25 26]);
RLCc{9} = int16([8 14 19 20 21 22 23 24 25 26]);
RLCc{10} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{11} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{12} = int16([9 14 19 20 21 22 23 24 25 26]);
RLCc{13} = int16([9 18 19 20 21 22 23 24 25 26]);
RLCc{14} = int16([11 18 19 20 21 22 23 24 25 26]);
RLCc{15} = int16([13 18 19 20 21 22 23 24 25 26]);
RLCc{16} = int16([11 17 18 19 20 21 22 23 24 25 26]);

switch C
    case 'Y'
        RLC = RLCy;
    case 'C'
        RLC = RLCc;
end
x1 = ZigZag(x);

k=2; Count = 0; Bits = double(0);
while k<= 64
    if x1(k) == 0
        Count = Count + 1;
        if k == 64
            Bits = Bits + double(RLC{1}(1));
        break;
        end
    else
        if Count == 0
            RL = Count;
            Level = round(log2(abs(x1(k)))+0.5);
            Bits = Bits + double(RLC{RL+1}(Level+1));
        elseif Count >=1 && Count <=15
            RL = Count;
            Level = round(log2(abs(x1(k)))+0.5);
            Bits = Bits + double(RLC{RL+1}(Level));
            Count = 0;
        else
            Bits = Bits + double(RLC{16}(1));
            Count = Count - 16;
        end
    end
    k = k+1;
end
end

function y = ZigZag(x)
  % Returneaza adresele scanate in zig-zag ale coef 8x8 dintr-o matrice
  y(1) = x(1, 1); y(2) = x(1, 2); y(3) = x(2, 1); y(4) = x(3, 1);
  y(5) = x(2, 2); y(6) = x(1, 3); y(7) = x(1, 4); y(8) = x(2, 3);
  y(9) = x(3, 2); y(10) = x(4, 1); y(11) = x(5, 1); y(12) = x(4, 2);
  y(13) = x(3, 3); y(14) = x(2, 4); y(15) = x(1, 5); y(16) = x(1, 6);
  y(17) = x(2, 5); y(18) = x(3, 4); y(19) = x(4, 3); y(20) = x(5, 2);
  y(21) = x(6, 1); y(22) = x(7, 1); y(23) = x(6, 2); y(24) = x(5, 3);
  y(25) = x(4, 4); y(26) = x(3, 5); y(27) = x(2, 6); y(28) = x(1, 7);
  y(29) = x(1, 8); y(30) = x(2, 7); y(31) = x(3, 6); y(32) = x(4, 5);
  y(33) = x(5, 4); y(34) = x(6, 3); y(35) = x(7, 2); y(36) = x(8, 1);
  y(37) = x(8, 2); y(38) = x(7, 3); y(39) = x(6, 4); y(40) = x(5, 5);
  y(41) = x(4, 6); y(42) = x(3, 7); y(43) = x(2, 8); y(44) = x(3, 8);
  y(45) = x(4, 7); y(46) = x(5, 6); y(47) = x(6, 5); y(48) = x(7, 4);
  y(49) = x(8, 3); y(50) = x(8, 4); y(51) = x(7, 5); y(52) = x(6, 6);
  y(53) = x(5, 7); y(54) = x(4, 8); y(55) = x(5, 8); y(56) = x(6, 7);
  y(57) = x(7, 6); y(58) = x(8, 5); y(59) = x(8, 6); y(60) = x(7, 7);
  y(61) = x(6, 8); y(62) = x(7, 8); y(63) = x(8, 7); y(64) = x(8, 8);
end

---

# Laborator 4 SCCV

### Aplicatia 1 (L4 SCCV)

clear all, close all, clc;

A = [2 2; -1 1];
[U S V] = svd(A);

### Aplicatia 2 (L4 SCCV)

clear all, close all, clc;

imfinfo('lena_256x256_BW.bmp')
[A, map] = imread ('lena_256x256_BW.bmp');
Ad = im2double (A);

figure('Name', 'Imaginea originala'), imshow(A);
[U, S, V] = svd(Ad);

C = zeros(size (Ad));
k = 5;
for j = 1:k
    C = C + U(:, j) * S(j, j) * V(:, j)';
end

X = uint8(255*C);

figure('Name', 'Imaginea reconstruita'), imshow (X);

filename = fullfile(pwd, ['LenaComprimata', num2str(k), '.bmp']);
imwrite(X, filename)

### Aplicatia 3 (L4 SCCV)
(Melania - testat)

clear all, close all, clc;

%Se citeste imaginea color Lena_color.png
figure('Name', 'Imaginea originala RGB')
L=imread('lena_color_512x512.bmp');
imshow(L)

%Se descompune imaginea in cele trei plane de culoare corespunzătoare
L1=L(:,:,1);
L2=L(:,:,2);
L3=L(:,:,3);
%Se convertesc valorile pixelilor din int in dubla precizie
I1=im2double(L1);
I2=im2double(L2);
I3=im2double(L3);

%Se afiseaza cele trei plane de culoare
figure('Name','Componenta de roşu')
imshow(I1)
figure('Name','Componenta de verde')
imshow(I2)
figure('Name','Componenta de albastru')
imshow(I3)

%Se calculează descompunerea SVD pentru cele trei matrice
[u1,s1,v1]=svd(I1);
[u2,s2,v2]=svd(I2);
[u3,s3,v3]=svd(I3);

%Se creaza trei matrice Cl, C2 si C3 ce se umplu initial cu zerouri
C1=zeros(size(I1));
C2=zeros(size(I2));
C3=zeros(size(I3));

%Setarea numărului de valori singulare nenule considerate pentru
%descompunere
k=100;

%Se reconstruiesc imaginile pe baza valorilor singulare
for j=1:k
    C1=C1+u1(:,j)*s1(j,j)*v1(:,j).';
    C2=C2+u2(:,j)*s2(j,j)*v2(:,j).';
    C3=C3+u3(:,j)*s3(j,j)*v3(:,j).';
end

%Afisarea imaginilor comprimate pe fiecare componenta de culoare
figure('Name','Componenta de roşu reconstruita')
imshow(C1)
figure('Name','Componenta de verde reconstruita')
imshow(C2)
figure('Name','Componenta de albastru reconstruita')
imshow(C3)

%Se convertesc valorile matricelor la intregi, fara semn reprezentaţi pe
%8 bits
R1=im2uint8(C1);
R2=im2uint8(C2);
R3=im2uint8(C3);

%Se creeaza o imagine Q ce are cate o componenta de culoare reprezentata de
%matricele Rl, R2, R3
Q(:,:,1)=R1;
Q(:,:,2)=R2;
Q(:,:,3)=R3;

%Afisarea imaginii finale comprimata
figure('Name', ' Imaginea color reconstruita')
imshow(Q,[])

---

# Laborator 5 SCCV

Melania
Va rog verificati! Am adaugat un "end" la sfarsit!

> `./main.m`

clear
close all
clc

A = imread('rino1.bmp'); % frame #1 - cadrul de referinţa
B = imread('rino2.bmp'); % frame #2 cadrul curent

N = 8; % Blocuri de dimensiune N x N pixeli
W = 16; % Dimensiunea ferestrei de cautare W x W pixeli

[x,y,ep,BRec] = MC_CautareCompleta(A,B,N,W);

C = double(B) - double(A);
figure, imshow(ep,[]), title('Predictia cadrului cu compensare de miscare')
figure, imshow(C, []), title('Predictia cadrului fara compensare de miscare')

> `./MC_CautareCompleta.m`

function[x, y, ep, BRec] = MC_CautareCompleta(A,B,N,W)
    % funcţia calculează vectorii de mişcare pentru fiecare NxN bloc din
    % imagine folosind tehnica de cautare completa
    % x - componenta orizontala a vectorului de mişcare
    % y - componenta verticala a vectorului de mişcare
    % ep - eroarea de predictie
    % BRec - imaginea reconstruita pe baza vectorilor de mişcare

    [Height,Width] = size (A);
    % se completează cadrul de referinţei la stanga, dreapta, sus si jos cu
    % pixeli pentru a ne asigura ca blocul de analiza se gaseste Întotdeauna in
    % centrul ferestrei de cautare, chiar si pentru blocurile aflate pe margine
    A1 = double(padarray(A,[(W-N)/2 (W-N)/2], 'replicate' )); 
    B1 = double(B);

    x = (zeros(Height/N,Width/N));
    y = (zeros(Height/N,Width/N));
    Xpoz = (zeros(Height/N,Width/N));
    Ypoz = (zeros(Height/N,Width/N));

    figure, imshow(B), title('Vectorii de miscare')
    hold on

    t1 = cputime; % start CPU time
    for r = 1:N:Height
        rblk = floor(r/N) + 1;
        for c = 1:N:Width
            cblk = floor(c/N) + 1;
            D = 1.0e+10;% distanta maxima intre doua blocuri
            for u = -(W-N)/2:(W-N)/2
                for v = -(W-N)/2:(W-N)/2
                    d = B1(r:r+N-1,c:c+N-1)-A1(r+u+(W-N)/2:r+u+(W-N)/2+N-1,c+v+(W-N)/2:c+v+(W-N)/2+N-1);
                    d = sum(abs(d(:)));
                    if d < D
                        D = d;
                        x(rblk, cblk) = v;
                        y(rblk, cblk);
                    end
                end
            end
            Xpoz(rblk, cblk) = c + N/2;
            Ypoz(rblk, cblk) = r + N/2;
        end
    end
    t2 = cputime; % end CPU time

    quiver(Xpoz, Ypoz, x, y, 'r','LineWidth', 1);
    hold off

    % Reconstruieşte cadrul curent utilizând cadrul de referinţa si ve%orii de
    % mişcare
    BRec = double(zeros(Height,Width)); % cadrul reconstruit
    ep = double(zeros(Height,Width)); % eroarea de predictie

    for r = 1:N:Height
        rblk = floor(r/N) + 1;
        for c=1:N:Width
            cblk = floor(c/N) + 1;
            x1 = x(rblk, cblk);
            y1 = y(rblk, cblk);
            BRec(r:r+N-1,c:c+N-1) = A1(r+y1+(W-N)/2:r+y1+(W-N)/2+N-1,c+x1+(W-N)/2:c+x1+(W-N)/2+N-1);
            ep(r:r+N-1,c:c+N-1) = B1(r:r+N-1,c:c+N-1) - BRec(r:r+N-1,c:c+N-1);
        end
    end
end