Eliminate block process.

Move parallelization into run script.

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

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

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

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