% Run Detection
thresholdC1 = 35.0;
thresholdC2 = 25.0;
thresholdU = 0.99;

for frameCnt = 1:size(ag,3)
    fprintf('Frame # %d\n', frameCnt);
    Nt = 0;
    frameWindow = ag(:,:,frameCnt:frameCnt+4);
    
    img = ag(:,:,frameCnt);
    for zy = 1:30
        for zx = 1:40
            zz = ((zy -1) * 40) + zx;
            yoffset = (zy - 1) * 8;
            xoffset = (zx - 1) * 8;

            i3Drow = 1;
            i3Dpix = 1;
            for i = 1:3
                for j = 1:3
                    data = frameWindow(:,:,i + (j - 1));
                    for py = 1:8
                        for px = 1:8
                            matrix3d(i3Drow, i3Dpix) = data(yoffset + py, xoffset + px);
                            i3Dpix = i3Dpix + 1;
                        end
                    end
                end
                i3Drow = i3Drow + 1;
            end

            for k = 1:3
                for l = 1:3
                    projected(k, l) = zeros(1, 1);
                    coloffset = BLOCK3D - l + 1;
                    for m = 1:BLOCK3D
                        projected(k, l) = projected(k, l) + (matrix3d(k, m) * projPCA{zz}(m, coloffset));
                    end
                end
            end

            [V D] = eig(cov(projected));
            ev = max(D(:));
            
            if((ev - blockmean(zz)) > (thresholdC1 * blockstd(zz)))
                Nt = Nt + 1;
                fprintf('Motion detected at region %d\n', zz);
            else
                stdupdate = sqrt(power((ev - blockmean(zz)), 2));
                blockmean(zz) = (thresholdU * blockmean(zz)) + (1.0 - thresholdU) * ev;
                blockstd(zz) = (thresholdU * blockstd(zz)) + (1.0 - thresholdU) * stdupdate;
            end
        end
    end
    subplot(1,2,1); imshow(ag(:,:,frameCnt)); title('Control/Original');
    subplot(1,2,2); imshow(img); title(sprintf('Number of Regions with Motion %d', Nt));
    pause(0.000005);
end