diff --git a/pcad.m b/pcad.m
new file mode 100644
index 0000000..35d5c7f
--- /dev/null
+++ b/pcad.m
@@ -0,0 +1,13 @@
+function im_result = pcad(img)
+%PCAD Perform the PCAD algorithm on img.
+
+height = size(img, 1);
+width = size(img, 2);
+
+cb = @(block) pcad_block(block, 50);
+r = blockproc(img, [15 15], cb, 'PadPartialBlocks', true, 'PadMethod', 'symmetric');
+    
+% MAX
+im_result = r(1:height, 1:width);
+
+end
diff --git a/pcad_block.m b/pcad_block.m
new file mode 100644
index 0000000..36f8f00
--- /dev/null
+++ b/pcad_block.m
@@ -0,0 +1,27 @@
+function px = pcad_block(block_struct, d)
+%PCAD_BLOCK Perform the PCAD algorithm on img
+
+% square dimension
+block_size = size(block_struct.data);
+channels = size(block_struct.data, 3);
+
+% PCA
+vals = reshape(block_struct.data, [], channels);
+[~, p_sc, p_la] = pca(vals, 'Algorithm', 'svd', 'Centered', false);
+
+% resuable values
+pad = nan(d, 1);
+p_nm = size(p_sc, 2);
+ret = zeros(size(p_sc));
+
+for c = 1:p_nm
+    [s, s_i] = sort(p_sc(:, c));
+    
+    ret(s_i, c) = max(s - [pad; s(1:end-d)], [s(1 + d:end); pad] - s);
+end
+
+% Euclidian distance
+% restore shape
+px = reshape(sqrt(sum(ret .^ 2, 2)), block_size(1), block_size(2));
+
+end
diff --git a/run_script.m b/run_script.m
index 0fb7470..b3eef40 100644
--- a/run_script.m
+++ b/run_script.m
@@ -97,6 +97,13 @@ function run_script(scene_file)
         save(fname, 'img_mwpcag');
     end
 
+    % PCAD
+    fname = sprintf('output/%s-%d-pcad.mat', scene_file, j);
+    if ~exist(fname, 'file')
+        img_pcad = pcad(img);
+        save(fname, 'img_pcad');
+    end
+
     % KNNA
     fname = sprintf('output/%s-%d-knna.mat', scene_file, j);
     if ~exist(fname, 'file')