diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..2fe6a26
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.tiff
+*.mat
diff --git a/colorspaces.txt b/colorspaces.txt
new file mode 100644
index 0000000..40cca91
--- /dev/null
+++ b/colorspaces.txt
@@ -0,0 +1,10 @@
+Color space transformations:
+
+function img_ss = select_channel(img, channels)
+img_ss = img(:, :, channels)
+end
+
+
+rgb2ycbcr
+@(img) rgb2ycbcr;
+im_rgb2lab
diff --git a/dwest.m b/dwest.m
new file mode 100644
index 0000000..896b634
--- /dev/null
+++ b/dwest.m
@@ -0,0 +1,7 @@
+function im_result = dwest(img)
+%DWEST Perform the DWEST algorithm on img.
+
+  
+im_result = blockproc(img, [1 1], @dwest_block, 'BorderSize', [5 5], 'TrimBorder', false);
+end
+
diff --git a/dwest_block.m b/dwest_block.m
new file mode 100644
index 0000000..77b1362
--- /dev/null
+++ b/dwest_block.m
@@ -0,0 +1,37 @@
+function px = dwest_block(block_struct)
+%DWEST Perform the DWEST algorithm on img.
+
+% square odd dimension
+block_size = size(block_struct.data, 1);
+channels = size(block_struct.data, 3);
+mid_point = (block_size - 1) / 2;
+mask = false([block_size block_size channels]);
+
+px = 0;
+for iw_size = 1:2:(block_size - 4)
+    % update mask
+    d = (iw_size - 1) / 2;
+    mask(mid_point - d:mid_point + d, mid_point - d:mid_point + d, :) = true([iw_size iw_size channels]);
+    
+    % extract windows
+    iw = reshape(block_struct.data(mask), [], channels);
+    ow = reshape(block_struct.data(~mask), [], channels);
+    
+    % calculate means
+    iw_mean = mean(iw, 1);
+    ow_mean = mean(ow, 1);
+    
+    % centered data
+    iw_centered = bsxfun(@minus, iw, iw_mean);
+    ow_centered = bsxfun(@minus, ow, ow_mean);
+    
+    % calculate covariance
+    iw_cov = (iw_centered' * iw_centered) / (size(ow_centered, 1) - 1);
+    ow_cov = (ow_centered' * ow_centered) / (size(ow_centered, 1) - 1);
+    
+    % eigenvalues difference in covariance
+    [e_vec, e_val] = eig(iw_cov - ow_cov);
+    
+    % eigen vectors associated with positive value times mean
+    px = max(abs(sum((ow_mean - iw_mean) * e_vec(:, diag(e_val) > 0))), px);
+end
diff --git a/im_rgb2lab.m b/im_rgb2lab.m
new file mode 100644
index 0000000..a15758b
--- /dev/null
+++ b/im_rgb2lab.m
@@ -0,0 +1,6 @@
+function im_lab = im_rgb2lab(im_rgb )
+%IM_RGB2LAB Convert RGB to L*a*b colorspace
+cf = makecform('srgb2lab');
+im_lab = applycform(im_rgb, cf);
+end
+
diff --git a/im_rgb2upvpl.m b/im_rgb2upvpl.m
new file mode 100644
index 0000000..2c3b94f
--- /dev/null
+++ b/im_rgb2upvpl.m
@@ -0,0 +1,7 @@
+function im_upvpl = im_rgb2upvpl(im_rgb )
+%IM_RGB2UPVPL Convert RGB to u'v'L colorspace
+cf = makecform('srgb2xyz');
+cf2 = makecform('xyz2upvpl');
+im_upvpl = applycform(applycform(im_rgb, cf), cf2);
+end
+
diff --git a/im_rgb2uvl.m b/im_rgb2uvl.m
new file mode 100644
index 0000000..68605fe
--- /dev/null
+++ b/im_rgb2uvl.m
@@ -0,0 +1,7 @@
+function im_uvl = im_rgb2uvl(im_rgb )
+%IM_RGB2LAB Convert RGB to uvL colorspace
+cf = makecform('srgb2xyz');
+cf2 = makecform('xyz2uvl');
+im_uvl = applycform(applycform(im_rgb, cf), cf2);
+end
+
diff --git a/im_rgb2xyy.m b/im_rgb2xyy.m
new file mode 100644
index 0000000..238e0dc
--- /dev/null
+++ b/im_rgb2xyy.m
@@ -0,0 +1,7 @@
+function im_xyy = im_rgb2xyy(im_rgb )
+%IM_RGB2XYY Convert RGB to xyY colorspace
+cf = makecform('srgb2xyz');
+cf2 = makecform('xyz2xyl');
+im_xyy = applycform(applycform(im_rgb, cf), cf2);
+end
+
diff --git a/im_rgb2xyz.m b/im_rgb2xyz.m
new file mode 100644
index 0000000..1d01625
--- /dev/null
+++ b/im_rgb2xyz.m
@@ -0,0 +1,6 @@
+function im_xyz = im_rgb2xyz(im_rgb)
+%IM_RGB2XYZ Convert RGB to XYZ colorspace
+cf = makecform('srgb2xyz');
+im_xyz = applycform(im_rgb, cf);
+end
+
diff --git a/nswtd.m b/nswtd.m
new file mode 100644
index 0000000..4920e0d
--- /dev/null
+++ b/nswtd.m
@@ -0,0 +1,6 @@
+function im_result = nswtd(img)
+%NSWTD Perform the NSWTD algorithm on img.
+  
+im_result = blockproc(img, [1 1], @nswtd_block, 'BorderSize', [5 5], 'TrimBorder', false);
+end
+
diff --git a/nswtd_block.m b/nswtd_block.m
new file mode 100644
index 0000000..b4b24e0
--- /dev/null
+++ b/nswtd_block.m
@@ -0,0 +1,30 @@
+function px = nswtd_block(block_struct)
+%NSWTD Perform the nested window anomaly detection algorithm.
+% Called for a single block of the iamge.
+
+% square odd dimension
+block_size = size(block_struct.data, 1);
+channels = size(block_struct.data, 3);
+mid_point = (block_size - 1) / 2;
+mask = false([block_size block_size channels]);
+
+px = 0;
+for iw_size = 1:2:(block_size - 4)
+    % update mask
+    d = (iw_size - 1) / 2;
+    mask(mid_point - d:mid_point + d, mid_point - d:mid_point + d, :) = true([iw_size iw_size channels]);
+    
+    % extract windows
+    iw = reshape(block_struct.data(mask), [], channels);
+    ow = reshape(block_struct.data(~mask), [], channels);
+    
+    % calculate means
+    iw_mean = mean(iw, 1);
+    ow_mean = mean(ow, 1);
+    
+    % OPD
+    d = opd(iw_mean, ow_mean);
+    
+    % eigen vectors associated with positive value times mean
+    px = max(d, px);
+end
diff --git a/opd.m b/opd.m
new file mode 100644
index 0000000..e287dce
--- /dev/null
+++ b/opd.m
@@ -0,0 +1,12 @@
+function d = opd(vec_a, vec_b)
+%OPD Orthogonal projection divergence
+
+d = (vec_a * (eye(numel(vec_a)) - vec_a' / (vec_a * vec_a') * vec_a) * vec_a') + (vec_b * (eye(numel(vec_b)) - vec_b' / (vec_b * vec_b') * vec_b) * vec_b');
+if 0 < d
+    d = sqrt(d);
+else
+    d = nan;
+end
+
+end
+
diff --git a/pcag.m b/pcag.m
new file mode 100644
index 0000000..e3eadaa
--- /dev/null
+++ b/pcag.m
@@ -0,0 +1,11 @@
+function im_result = pcag(img, sz)
+%PCAG Perform the PCAG algorithm on img.
+  
+if nargin < 2
+    sz = 5;
+end
+
+im_result = blockproc(img, [sz sz], @pcag_block);
+
+end
+
diff --git a/pcag_block.m b/pcag_block.m
new file mode 100644
index 0000000..750662d
--- /dev/null
+++ b/pcag_block.m
@@ -0,0 +1,21 @@
+function px = pcag_block(block_struct)
+%PCAG Perform the PCAG algorithm on img.
+
+% square odd dimension
+block_size = size(block_struct.data);
+channels = size(block_struct.data, 3);
+
+vals = reshape(block_struct.data, [], channels);
+[~, p_sc, p_la] = pca(vals);
+
+n = 0;
+p_nm = size(p_sc, 2);
+rm = ((block_size(1) * block_size(2)) - ((block_size(1) - 2) * (block_size(2) - 2))) / 2;
+for c = 1:p_nm
+    s = sort(p_sc(:, c));
+    n = n + max(diff(s(rm:end - rm))) * p_la(c);
+end
+n = n / sum(p_la);
+
+px = n * ones(block_size(1:2));
+end
diff --git a/sq_mask.m b/sq_mask.m
new file mode 100644
index 0000000..5e37180
--- /dev/null
+++ b/sq_mask.m
@@ -0,0 +1,18 @@
+function mask = sq_mask(sz, channels, m_size, h_size)
+%SQ_MASK Generate a square mask.
+% sz must be an odd number represnting the mask size
+mid_point = (sz - 1) / 2;
+if m_size == sz
+    mask = true([sz sz channels]);
+else
+    mask = false([sz sz channels]);
+    d = (m_size - 1) / 2;
+    mask(mid_point - d:mid_point + d, mid_point - d:mid_point + d, :) = true;
+end
+% hollow
+if nargin == 4
+    d = (h_size - 1) / 2;
+    mask(mid_point - d:mid_point + d, mid_point - d:mid_point + d, :) = false;
+end
+end
+    
\ No newline at end of file