From c03359aa8e91c7092e524ee606baf72885517885 Mon Sep 17 00:00:00 2001
From: Nathan Perkins <nathan@maxmo.net>
Date: Sun, 5 Apr 2015 17:01:25 -0400
Subject: [PATCH] Eliminate block process. Move parallelization into run
 script.

---
 dwest.m               | 62 +++++++++++++++++++++++++++++++++-
 dwest_block.m         | 34 -------------------
 mwnswtd.m             | 58 ++++++++++++++++++++++++++++++--
 mwnswtd_block.m       | 41 -----------------------
 nswtd.m               | 48 ++++++++++++++++++++++++--
 nswtd_block.m         | 30 -----------------
 pad_image_symmetric.m | 33 ++++++++++++++++++
 run.m                 | 71 ++-------------------------------------
 run_script.m          | 78 +++++++++++++++++++++++++++++++++++++++++++
 9 files changed, 276 insertions(+), 179 deletions(-)
 delete mode 100644 dwest_block.m
 delete mode 100644 mwnswtd_block.m
 delete mode 100644 nswtd_block.m
 create mode 100644 pad_image_symmetric.m
 create mode 100644 run_script.m

diff --git a/dwest.m b/dwest.m
index 8bea7d2..c63811a 100644
--- a/dwest.m
+++ b/dwest.m
@@ -1,4 +1,64 @@
 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, 'PadMethod', 'symmetric', 'UseParallel', true);
+
+% largest window size (must be odd)
+w_max_size = 11;
+
+% image information
+[height, width, channels] = size(img);
+
+% return image
+im_result = zeros(height, width);
+
+% prebuild masks
+masks = {};
+masks_inv = {};
+for w_size = 1:2:(w_max_size - 2)
+	masks{end + 1} = reshape(sq_mask(w_max_size, 1, w_size), [], 1);
+	masks_inv{end + 1} = ~masks{end};
+end
+num_masks = length(masks);
+
+% pad image
+border = (w_max_size - 1) / 2;
+pad_img = pad_image_symmetric(img, border);
+for i = 1:width
+	for j = 1:height
+		d = nan;
+		
+		% reshape window
+		win = reshape(pad_img(j:j + w_max_size - 1, i:i + w_max_size - 1, :), [], channels);
+		
+		% for each mask
+		for k = 1:num_masks
+			iw = win(masks{k}, :);
+			ow = win(masks_inv{k}, :);
+    
+			% inner window
+            if 1 < k
+                iw_mean = mean(iw, 1);
+                iw_centered = bsxfun(@minus, iw, iw_mean);
+            else
+                iw_mean = iw;
+                iw_centered = zeros(size(iw));
+            end
+			
+            % outer window
+			ow_mean = mean(ow, 1);
+			ow_centered = bsxfun(@minus, ow, ow_mean);
+    
+			% difference in covariance
+            diff_cov = ((iw_centered' * iw_centered) - (ow_centered' * ow_centered)) / (size(ow_centered, 1) - 1);
+    
+			% eigenvalues difference in covariance
+			[e_vec, e_val] = eig(diff_cov);
+
+			% eigen vectors associated with positive value times mean
+			d = max(d, abs(sum((ow_mean - iw_mean) * e_vec(:, diag(e_val) > 0))));
+		end
+		
+		im_result(j, i) = d;
+	end
+end
+
 end
diff --git a/dwest_block.m b/dwest_block.m
deleted file mode 100644
index 4919a72..0000000
--- a/dwest_block.m
+++ /dev/null
@@ -1,34 +0,0 @@
-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);
-
-px = nan;
-for iw_size = 1:2:(block_size - 4)
-    % update mask
-    mask = sq_mask(block_size, channels, iw_size);
-    
-    % 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/mwnswtd.m b/mwnswtd.m
index adada30..98f6698 100644
--- a/mwnswtd.m
+++ b/mwnswtd.m
@@ -1,6 +1,60 @@
 function im_result = mwnswtd(img)
 %MWNSWTD Perform the MW-NSWTD algorithm on img.
-  
-im_result = blockproc(img, [1 1], @mwnswtd_block, 'BorderSize', [5 5], 'TrimBorder', false, 'PadMethod', 'symmetric', 'UseParallel', true);
+
+% largest window size (must be odd)
+w_max_size = 11;
+
+% image information
+[height, width, channels] = size(img);
+
+% return image
+im_result = zeros(height, width);
+
+% create inner window (size 1)
+iw_size = 1;
+mask_iw = reshape(sq_mask(w_max_size, 1, iw_size), [], 1);
+
+% prebuild other masks
+masks = {};
+for w_size = (iw_size + 2):2:w_max_size
+	masks{end + 1} = reshape(sq_mask(w_max_size, 1, w_size, w_size - 2), [], 1);
+end
+num_masks = length(masks);
+
+% pad image
+border = (w_max_size - 1) / 2;
+pad_img = pad_image_symmetric(img, border);
+for i = 1:width
+	for j = 1:height
+		d = nan;
+		
+		% reshape window
+		win = reshape(pad_img(j:j + w_max_size - 1, i:i + w_max_size - 1, :), [], channels);
+		
+		% get inner window mean
+		iw_mean = mean(win(mask_iw, :), 1);
+        
+        % last outer window to pass into loop
+        ow_mean = mean(win(masks{1}, :), 1);
+		
+		% for each mask
+		for k = 2:num_masks
+			% last outer window is new middle window
+            mw_mean = ow_mean;
+            
+            % new outer window
+            ow_mean = mean(win(masks{k}, :), 1);
+	
+            % use OSP based on outer window to project both inner and middle window
+            p_outer = osp(ow_mean);
+            cd = (iw_mean * p_outer * iw_mean') + (mw_mean * p_outer * mw_mean');
+            if 0 <= cd
+                d = max(d, sqrt(cd));
+            end
+		end
+		
+		im_result(j, i) = d;
+	end
 end
 
+end
diff --git a/mwnswtd_block.m b/mwnswtd_block.m
deleted file mode 100644
index 512d09a..0000000
--- a/mwnswtd_block.m
+++ /dev/null
@@ -1,41 +0,0 @@
-function px = mwnswtd_block(block_struct)
-%MWNSWTD_BLOCK Perform the multi-window nested window anomaly detection
-%algorithm.
-%   Called for a single block fof the image.
-
-% square odd dimension
-block_size = size(block_struct.data, 1);
-channels = size(block_struct.data, 3);
-
-% create inner window (size 1)
-iw_size = 1;
-mask = sq_mask(block_size, channels, iw_size);
-iw_mean = mean(reshape(block_struct.data(mask), [], channels), 1);
-
-% mask to pass into loop
-mask = sq_mask(block_size, channels, iw_size + 2, iw_size);
-ow_mean = mean(reshape(block_struct.data(mask), [], channels), 1);
-
-px = nan;
-for w_size = (iw_size + 2):2:(block_size - 2)
-	% mask is already correct for middle window
-	mw_mean = ow_mean;
-	
-	% redraw mask for outerwindow
-	mask = sq_mask(block_size, channels, w_size + 2, w_size);
-	ow_mean = mean(reshape(block_struct.data(mask), [], channels), 1);
-	
-    % use OSP based on outer window to project both inner and middle window
-    p_outer = osp(ow_mean);
-    d = (iw_mean * p_outer * iw_mean') + (mw_mean * p_outer * mw_mean');
-    if 0 <= d
-        d = sqrt(d);
-    else
-        d = nan;
-    end
-    
-    % eigen vectors associated with positive value times mean
-    px = max(d, px);
-end
-
-end
diff --git a/nswtd.m b/nswtd.m
index 2902573..0614782 100644
--- a/nswtd.m
+++ b/nswtd.m
@@ -1,6 +1,50 @@
 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, 'PadMethod', 'symmetric', 'UseParallel', true);
+
+% largest window size (must be odd)
+w_max_size = 11;
+
+% image information
+[height, width, channels] = size(img);
+
+% return image
+im_result = zeros(height, width);
+
+% prebuild masks
+masks = {};
+masks_inv = {};
+for w_size = 1:2:(w_max_size - 2)
+	masks{end + 1} = reshape(sq_mask(w_max_size, 1, w_size), [], 1);
+	masks_inv{end + 1} = ~masks{end};
 end
+num_masks = length(masks);
 
+% pad image
+border = (w_max_size - 1) / 2;
+pad_img = pad_image_symmetric(img, border);
+for i = 1:width
+	for j = 1:height
+		d = nan;
+		
+		% reshape window
+		win = reshape(pad_img(j:j + w_max_size - 1, i:i + w_max_size - 1, :), [], channels);
+		
+		% for each mask
+		for k = 1:num_masks
+            % inner and outer window mean
+            if 1 < k
+                iw_mean = mean(win(masks{k}, :), 1);
+            else
+                iw_mean = win(masks{k}, :);
+            end
+			ow_mean = mean(win(masks_inv{k}, :), 1);
+            
+            % OPD
+            d = opd(iw_mean, ow_mean);
+		end
+		
+		im_result(j, i) = d;
+	end
+end
+
+end
diff --git a/nswtd_block.m b/nswtd_block.m
deleted file mode 100644
index ec3f40a..0000000
--- a/nswtd_block.m
+++ /dev/null
@@ -1,30 +0,0 @@
-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 = nan;
-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/pad_image_symmetric.m b/pad_image_symmetric.m
new file mode 100644
index 0000000..e2f7861
--- /dev/null
+++ b/pad_image_symmetric.m
@@ -0,0 +1,33 @@
+function img_padded = pad_image_symmetric(img, padding)
+%PAD_IMAGE_SYMMETRIC Pad img by padding semmetrically filling in new pixels
+
+% constant padding
+if all(size(padding) == 1)
+    padding = [padding padding padding padding];
+end
+
+% expand padding
+pad_top = padding(1);
+pad_right = padding(1);
+pad_bottom = padding(3);
+pad_left = padding(4);
+
+% get current dimensions
+[height, width, channels] = size(img);
+
+img_padded = zeros(height + pad_top + pad_bottom, width + pad_left + pad_right, channels);
+img_padded(pad_top + 1:pad_top + height, pad_left + 1:pad_left + width, :) = img;
+
+% fill in sides
+img_padded(1:pad_top, pad_left + 1:pad_left + width, :) = img(pad_top:-1:1, :, :);
+img_padded(pad_top + height:end, pad_left + 1:pad_left + width, :) = img(end:-1:end - pad_bottom, :, :);
+img_padded(pad_top + 1:pad_top + height, 1:pad_left, :) = img(:, pad_left:-1:1, :);
+img_padded(pad_top + 1:pad_top + height, pad_left + width:end, :) = img(:, end:-1:end - pad_right, :);
+
+% fill in corners
+img_padded(1:pad_top, 1:pad_left, :) = img(pad_top:-1:1, pad_left:-1:1, :);
+img_padded(1:pad_top, pad_left + width:end, :) = img(pad_top:-1:1, end:-1:end - pad_right, :);
+img_padded(pad_top + height:end, 1:pad_left, :) = img(end:-1:end - pad_bottom, pad_left:-1:1, :);
+img_padded(pad_top + height:end, pad_left + width:end, :) = img(end:-1:end - pad_bottom, end:-1:end - pad_right, :);
+
+end
diff --git a/run.m b/run.m
index 768a5b0..da33aec 100644
--- a/run.m
+++ b/run.m
@@ -2,79 +2,12 @@
 scene_files = dir('scenes/*.jpg');
 scene_files = {scene_files.name};
 
-% color spaces
-color_spaces = {};
-color_spaces{end + 1} = @(img) img;
-color_spaces{end + 1} = @im_rgb2lab;
-color_spaces{end + 1} = @im_rgb2upvpl;
-color_spaces{end + 1} = @im_rgb2uvl;
-color_spaces{end + 1} = @im_rgb2xyy;
-color_spaces{end + 1} = @im_rgb2xyz;
-color_spaces{end + 1} = @(img) select_channel(im_rgb2lab(img), 2:3);
-color_spaces{end + 1} = @(img) select_channel(im_rgb2upvpl(img), 1:2);
-color_spaces{end + 1} = @(img) select_channel(im_rgb2uvl(img), 1:2);
-color_spaces{end + 1} = @(img) select_channel(im_rgb2xyz(img), [1 3]);
-color_spaces{end + 1} = @(img) select_channel(im_rgb2lab(img), [1 3]);
-
 % tuned for my laptop (3 workers)
 parpool('local', 3);
 
 % build
-for i = 1:length(scene_files)
-    % file name
-    scene_file = scene_files{i};
-    
-    % build scene and target
-    [scene, target] = build_scene(scene_file);
-    
-    % save
-    save(['output/' scene_file '.mat'], 'scene', 'target');
-    
-    fprintf('Scene %s...\n', scene_file);
-    
-    % for each color space
-    for j = 1:length(color_spaces)
-        % get callback
-        color_cb = color_spaces{j};
-        
-        % transform image
-        img = color_cb(scene);
-        
-        fprintf('Color %d...\n', j);
-        
-        % run four algorithms
-        fname = sprintf('output/%s-%d-dwest.mat', scene_file, j);
-        if exist(fname, 'file')
-            load(fname);
-        else
-            img_dwest = dwest(img);
-            save(fname, 'img_dwest');
-        end
-        
-        fname = sprintf('output/%s-%d-nswtd.mat', scene_file, j);
-        if exist(fname, 'file')
-            load(fname);
-        else
-            img_nswtd = nswtd(img);
-            save(fname, 'img_nswtd');
-        end
-        
-        fname = sprintf('output/%s-%d-mwnswtd.mat', scene_file, j);
-        if exist(fname, 'file')
-            load(fname);
-        else
-            img_mwnswtd = mwnswtd(img);
-            save(fname, 'img_mwnswtd');
-        end
-        
-        fname = sprintf('output/%s-%d-pcag.mat', scene_file, j);
-        if exist(fname, 'file')
-            load(fname);
-        else
-            img_pcag = mwnswtd(img);
-            save(fname, 'img_pcag');
-        end
-    end 
+parfor i = 1:length(scene_files)
+    run_script(scene_files{i});
 end
 
 delete(gcp);
\ No newline at end of file
diff --git a/run_script.m b/run_script.m
new file mode 100644
index 0000000..4f6a1ca
--- /dev/null
+++ b/run_script.m
@@ -0,0 +1,78 @@
+function run_script(scene_file)
+%RUN_SCRIPT
+
+% color spaces
+color_spaces = {};
+color_spaces{end + 1} = @(img) img;
+color_spaces{end + 1} = @im_rgb2lab;
+color_spaces{end + 1} = @im_rgb2upvpl;
+color_spaces{end + 1} = @im_rgb2uvl;
+color_spaces{end + 1} = @im_rgb2xyy;
+color_spaces{end + 1} = @im_rgb2xyz;
+color_spaces{end + 1} = @(img) select_channel(im_rgb2lab(img), 2:3);
+color_spaces{end + 1} = @(img) select_channel(im_rgb2upvpl(img), 1:2);
+color_spaces{end + 1} = @(img) select_channel(im_rgb2uvl(img), 1:2);
+color_spaces{end + 1} = @(img) select_channel(im_rgb2xyz(img), [1 3]);
+color_spaces{end + 1} = @(img) select_channel(im_rgb2lab(img), [1 3]);
+color_spaces{end + 1} = @(img) transform_channel(im_rgb2lab(img), 1, @log);
+
+% load file
+fname = sprintf('output/%s.mat', scene_file);
+if exist(fname, 'file')
+    load(fname);
+else
+    % build scene and target
+    [scene, target] = build_scene(scene_file);
+
+    % save
+    save(['output/' scene_file '.mat'], 'scene', 'target');
+end
+
+fprintf('Scene %s...\n', scene_file);
+
+% for each color space
+for j = 1:length(color_spaces)
+    % get callback
+    color_cb = color_spaces{j};
+
+    % transform image
+    img = color_cb(scene);
+
+    fprintf('Color %d...\n', j);
+
+    % run four algorithms
+    fname = sprintf('output/%s-%d-dwest.mat', scene_file, j);
+    if exist(fname, 'file')
+        load(fname);
+    else
+        img_dwest = dwest(img);
+        save(fname, 'img_dwest');
+    end
+
+    fname = sprintf('output/%s-%d-nswtd.mat', scene_file, j);
+    if exist(fname, 'file')
+        load(fname);
+    else
+        img_nswtd = nswtd(img);
+        save(fname, 'img_nswtd');
+    end
+
+    fname = sprintf('output/%s-%d-mwnswtd.mat', scene_file, j);
+    if exist(fname, 'file')
+        load(fname);
+    else
+        img_mwnswtd = mwnswtd(img);
+        save(fname, 'img_mwnswtd');
+    end
+
+    fname = sprintf('output/%s-%d-pcag.mat', scene_file, j);
+    if exist(fname, 'file')
+        load(fname);
+    else
+        img_pcag = mwnswtd(img);
+        save(fname, 'img_pcag');
+    end
+end 
+
+end
+