Initial commit

LICENSE

@@ -0,0 +1,17 @@
+Copyright 2018, Yagiz Aksoy. All rights reserved.
+
+This software is for academic use only. A redistribution of this 
+software, with or  without modifications, has to be for academic 
+use only, while giving the appropriate credit to the original 
+authors of the software. The methods implemented as a part of 
+this software may be covered under patents or patent applications.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHOR ''AS IS'' AND ANY EXPRESS OR IMPLIED
+WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

SemanticSoftSegmentation.m

@@ -0,0 +1,49 @@
+
+% Semantic Soft Segmentation
+% This function implements the soft segmentation approach described in
+% Yagiz Aksoy, Tae-Hyun Oh, Sylvain Paris, Marc Pollefeys, Wojciech Matusik
+% "Semantic Soft Segmentation", ACM TOG (Proc. SIGGRAPH) 2018
+
+function [softSegments, initSoftSegments, Laplacian, affinities, features, superpixels, eigenvectors, eigenvalues] = SemanticSoftSegmentation(image, features)
+
+    disp('Semantic Soft Segmentation')
+    % Prepare the inputs and superpixels
+    image = im2double(image);
+    if size(features, 3) > 3 % If the features are raw, hyperdimensional, preprocess them
+        features = preprocessFeatures(features, image);
+    else
+        features = im2double(features);
+    end
+    superpixels = Superpixels(image);
+    [h, w, ~] = size(image);
+
+    disp('     Computing affinities')
+    % Compute the affinities and the Laplacian
+    affinities{1} = mattingAffinity(image);
+    affinities{2} = superpixels.neighborAffinities(features); % semantic affinity
+    affinities{3} = superpixels.nearbyAffinities(image); % non-local color affinity
+    Laplacian = affinityMatrixToLaplacian(affinities{1} + 0.01 * affinities{2} + 0.01 * affinities{3}); % Equation 6
+
+    disp('     Computing eigenvectors')
+    % Compute the eigendecomposition
+    eigCnt = 100; % We use 100 eigenvectors in the optimization
+    [eigenvectors, eigenvalues] = eigs(Laplacian, eigCnt, 'SM');
+
+    disp('     Initial optimization')
+    % Compute initial soft segments
+    initialSegmCnt = 40;
+    sparsityParam = 0.8;
+    iterCnt = 40;
+    % feeding features to the function below triggers semantic intialization
+    initSoftSegments = softSegmentsFromEigs(eigenvectors, eigenvalues, Laplacian, ...
+                                            h, w, features, initialSegmCnt, iterCnt, sparsityParam, [], []);
+
+    % Group segments w.r.t. their mean semantic feature vectors
+    groupedSegments = groupSegments(initSoftSegments, features);
+
+    disp('     Final optimization')
+    % Do the final sparsification
+    softSegments = sparsifySegments(groupedSegments, Laplacian, imageGradient(image, false, 6));
+
+    disp('     Done.')
+end

SpectralMatting.m

@@ -0,0 +1,31 @@
+
+% Spectral Matting
+% This function implements the soft segmentation approach described in
+% Anat Levin, Dani Lischinski, Yair Weiss, "Spectral Matting", IEEE TPAMI, 2008.
+% The parameters here are set to their default values as reported by Levin et al.
+
+function [softSegments, Laplacian, eigenvectors, eigenvalues] = SpectralMatting(image)
+
+    disp('Spectral Matting')
+    image = im2double(image);
+    [h, w, ~] = size(image);
+
+    disp('     Computing affinities')
+    % Compute the matting Laplacian
+    Laplacian = affinityMatrixToLaplacian(mattingAffinity(image));
+
+    disp('     Computing eigenvectors')
+    % Compute the eigendecomposition
+    eigCnt = 50;
+    [eigenvectors, eigenvalues] = eigs(Laplacian, eigCnt, 'SM');
+
+    disp('     Optimization')
+    % Compute the soft segments = matting components
+    initialSegmCnt = 20;
+    sparsityParam = 0.8;
+    iterCnt = 20;
+    softSegments = softSegmentsFromEigs(eigenvectors, eigenvalues, Laplacian, ...
+                                            h, w, [], initialSegmCnt, iterCnt, sparsityParam, [], []);
+
+    disp('     Done.')
+end

Superpixels.m

@@ -0,0 +1,125 @@
+
+% This class is for computing superpixel-based affinities described in the paper.
+% This class requires the image graphs methods by Steve Eddins. Find it here:
+% http://www.mathworks.com/matlabcentral/fileexchange/53614-image-graphs
+
+classdef Superpixels
+    properties 
+        labels
+        spcount
+        neigh
+        centroids
+    end
+    methods
+        function obj = Superpixels(im, spcnt)
+            if ~exist('spcnt', 'var') || isempty(spcnt)
+                spcnt = 2500;
+            end
+            [L, N] = superpixels(im, spcnt, 'Compactness', 1e-20);
+            obj.labels = L;
+            obj.spcount = N;
+            % Find neighboring superpixels
+            g = adjacentRegionsGraph(L);
+            obj.neigh = g.Edges.Labels;
+            % Find centroids
+            s = regionprops(L, 'centroid');
+            cent = cat(1, s.Centroid);
+            obj.centroids = round(cent(:, 2:-1:1));
+            [h, w, ~] = size(im);
+            obj.centroids(:, 3) = sub2ind([h, w], obj.centroids(:, 1), obj.centroids(:, 2));
+        end
+
+        function regmeans = computeRegionMeans(obj, image)
+            [h, w, c] = size(image);
+            image = reshape(image, [h*w, c]);
+            regmeans = zeros(obj.spcount, c);
+            idx = label2idx(obj.labels);
+            for i = 1 : length(idx)
+                regmeans(i, :) = mean(image(idx{i}, :), 1);
+            end
+        end
+
+        % This is for the semantic affinity, generates affinities in [-1, 1]
+        function W = neighborAffinities(obj, features, erfSteepness, erfCenter)
+            if ~exist('erfSteepness', 'var') || isempty(erfSteepness)
+                erfSteepness = 20;
+            end
+            if ~exist('erfCenter', 'var') || isempty(erfCenter)
+                erfCenter = 0.85;
+            end
+            [h, w, ~] = size(features);
+            N = h * w;
+            spMeans = obj.computeRegionMeans(features);
+            affs = zeros(size(obj.neigh, 1), 1);
+            inds1 = affs;
+            inds2 = affs;
+            for i = 1 : size(obj.neigh, 1)
+                ind1 = obj.neigh(i, 1);
+                ind2 = obj.neigh(i, 2);
+                affs(i) = sigmoidAff(spMeans(ind1, :), spMeans(ind2, :), erfSteepness, erfCenter);
+                inds1(i) = obj.centroids(ind1, 3);
+                inds2(i) = obj.centroids(ind2, 3);
+            end
+            W = sparse(inds1, inds2, affs, N, N);
+            W = W' + W;
+        end
+
+        % This is for the nonlocal color affinity, generates affinities in [0, 1]
+        function W = nearbyAffinities(obj, image, erfSteepness, erfCenter, proxThresh)
+            if ~exist('erfSteepness', 'var') || isempty(erfSteepness)
+                erfSteepness = 50;
+            end
+            if ~exist('erfCenter', 'var') || isempty(erfCenter)
+                erfCenter = 0.95;
+            end
+            if ~exist('proxThresh', 'var') || isempty(proxThresh)
+                proxThresh = 0.2;
+            end
+            [h, w, ~] = size(image);
+            N = h * w;
+            spMeans = obj.computeRegionMeans(image);
+            combinationCnt = obj.spcount;
+            combinationCnt = combinationCnt * (combinationCnt - 1) / 2;
+            affs = zeros(combinationCnt, 1);
+            inds1 = affs;
+            inds2 = affs;
+            cnt = 1;
+            cents = obj.centroids(:, 1:2);
+            cents(:,1) = cents(:,1) / h;
+            cents(:,2) = cents(:,2) / w;
+            for i = 1 : obj.spcount
+                for j = i + 1 : obj.spcount
+                    centdist = cents(i, 1:2) - cents(j, 1:2);
+                    centdist = sqrt(centdist * centdist');
+                    if centdist > proxThresh
+                        affs(cnt) = 0;
+                    else
+                        affs(cnt) = sigmoidAffPos(spMeans(i, :), spMeans(j, :), erfSteepness, erfCenter);
+                    end
+                    inds1(cnt) = obj.centroids(i, 3);
+                    inds2(cnt) = obj.centroids(j, 3);
+                    cnt = cnt + 1;
+                end
+            end
+            W = sparse(inds1, inds2, affs, N, N);
+            W = W' + W;
+        end
+
+        function vis = visualizeRegionMeans(obj, im)
+            vis = label2rgb(obj.labels, obj.computeRegionMeans(im));
+        end
+
+    end
+end
+
+function aff = sigmoidAff(feat1, feat2, steepness, center)
+    aff = abs(feat1 - feat2);
+    aff = 1 - sqrt(aff * aff');
+    aff = (erf(steepness * (aff - center)));
+end
+
+function aff = sigmoidAffPos(feat1, feat2, steepness, center)
+    aff = abs(feat1 - feat2);
+    aff = 1 - sqrt(aff * aff');
+    aff = (erf(steepness * (aff - center)) + 1) / 2;
+end

affinityMatrixToLaplacian.m

@@ -0,0 +1,16 @@
+
+function Lap = affinityMatrixToLaplacian(aff, normalize)
+    if ~exist('normalize', 'var') || isempty(normalize)
+        normalize = false ;
+    end
+    N = size(aff, 1);
+    if normalize
+        aa = sum(aff, 2);
+        D = spdiags(aa, 0 , N, N);
+        aa = sqrt(1./aa);
+        D12 = spdiags(aa, 0 , N, N);
+        Lap = D12 * (D - aff) * D12;
+    else
+        Lap = spdiags(sum(aff, 2), 0 , N, N) - aff;
+    end
+end

demo.m

@@ -0,0 +1,40 @@
+
+% Add the ImageGraph to path (in a folder named 'ImageGraphs'). Find it here:
+% http://www.mathworks.com/matlabcentral/fileexchange/53614-image-graphs
+addpath(fullfile(fileparts(mfilename('fullpath')), 'ImageGraphs'));
+
+%% Read the image and features from the sample file
+image = im2double(imread('docia.png'));
+features = image(:, size(image, 2) / 2 + 1 : end, :);
+image = image(:, 1 : size(image, 2) / 2, :);
+
+% The eigendecomposition uses a lot of memory and may render the computer
+% unresponsive, so better to test it first with a small image.
+image = imresize(image, 0.5);
+features = imresize(features, 0.5);
+
+%% Semantic soft segmentation
+% This function outputs many intermediate variables, if needed.
+% The results may vary a bit from run to run, as there are 2 stages that use 
+% k-means for intialization & grouping.
+sss = SemanticSoftSegmentation(image, features);
+
+% To use the features generated using our network implementation,
+% just feed them as the 'features' variable to the function. It will do
+% the prepocessing described in the paper and give the processed
+% features as an output.
+% If you are dealing with many images, storing the features after
+% preprocessing is recommended as raw hyperdimensional features
+% take a lot of space. Check the 'preprocessFeatures.m' file.
+
+% Visualize
+figure; imshow([image features visualizeSoftSegments(sss)]);
+title('Semantic soft segments');
+
+% There's also an implementation of Spectral Matting included
+sm = SpectralMatting(image);
+% You can group the soft segments from Spectral Matting using
+% semantic features, the way we presented our comparisons in the paper.
+sm_gr = groupSegments(sm, features);
+figure; imshow([image visualizeSoftSegments(sm) visualizeSoftSegments(sm_gr)]);
+title('Matting components');

groupSegments.m

@@ -0,0 +1,21 @@
+
+% A simple grouping of soft segments w.r.t. their semantic features
+% as described in Section 3.4.
+
+function groupedSegments = groupSegments(segments, features, segmCnt)
+    if ~exist('segmCnt', 'var') || isempty(segmCnt)
+        segmCnt = 5;
+    end
+    [h, w, cnt] = size(segments);
+    compFeatures = zeros(cnt, size(features, 3));
+    for i = 1 : cnt
+        cc = segments(:,:,i) .* features;
+        cc = sum(sum(cc, 1), 2) / sum(sum(segments(:,:,i), 1), 2);
+        compFeatures(i, :) = cc;
+    end
+    ids = kmeans(compFeatures, segmCnt);
+    groupedSegments = zeros(h, w, segmCnt);
+    for i = 1 : segmCnt
+        groupedSegments(:,:,i) = sum(segments(:,:,ids==i), 3);
+    end
+end

imageGradient.m

@@ -0,0 +1,67 @@
+
+% Yagiz Aksoy, 2016
+% Implements H. Farid, E.P. Simoncelli, Differentiation of Discrete Multidimensional Signals, TIP 2004
+% outColor switch makes the output 3-channel (color) or 1-channel, default true
+% There are 3 different filtSize options (1, 4 or 6), which determines the number of taps in the filters
+% and hence which neighborhood is used to approximate the derivatives
+
+function [gradientMagnitude, gradientOrientation, xGradient, yGradient] = imageGradient(image, outColor, filtSize)
+    % Set up variables legally
+    if ~exist('outColor', 'var') || isempty(outColor)
+        outColor = true;
+    end
+    if ~exist('filtSize', 'var') || isempty(filtSize)
+        filtSize = 1;
+    end
+    if filtSize <= 3
+        filtSize = 1;
+    elseif filtSize <= 5
+        filtSize = 4;
+    else
+        filtSize = 6;
+    end
+    image = im2double(image);
+
+    % Set up one-dimensional filters
+    switch filtSize
+        case 1
+            dk = [0.425287, -0.0000, -0.425287];
+            kk = [0.229879, 0.540242, 0.229879];
+        case 4
+            dk = [0.0032, 0.0350, 0.1190, 0.1458, -0.0000, -0.1458, -0.1190, -0.0350, -0.0032];
+            kk = [0.0009, 0.0151, 0.0890, 0.2349, 0.3201, 0.2349, 0.0890, 0.0151, 0.0009];
+        case 6
+            dk = [0.0001, 0.0019, 0.0142, 0.0509, 0.0963, 0.0878, 0.0000, -0.0878, -0.0963, -0.0509, -0.0142, -0.0019, -0.0001];
+            kk = [0.0000, 0.0007, 0.0071, 0.0374, 0.1126, 0.2119, 0.2605, 0.2119, 0.1126, 0.0374, 0.0071, 0.0007, 0.0000];
+    end
+
+    % Repeat-pad image
+    leftPad = image(:, 1, :);
+    rightPad = image(:, end, :);
+    image = [repmat(leftPad, [1 13 1]), image, repmat(rightPad, [1 13 1])];
+    upPad = image(1, :, :);
+    downPad = image(end, :, :);
+    image  = [repmat(upPad, [13 1 1]); image; repmat(downPad, [13 1 1])];
+
+    % Compute gradients
+    yGradient = zeros(size(image));
+    xGradient = zeros(size(image));
+    for i = 1 : size(image, 3)
+        yGradient(:,:,i) = conv2(dk, kk, image(:,:,i), 'same');
+        xGradient(:,:,i) = conv2(kk, dk, image(:,:,i), 'same');
+    end
+
+    % Remove padding
+    yGradient = yGradient(14 : end - 13, 14 : end - 13, :);
+    xGradient = xGradient(14 : end - 13, 14 : end - 13, :);
+
+    % Compute pixel-wise L2 norm if no color option is selected
+    if ~outColor
+        yGradient = sqrt(sum(yGradient .* yGradient, 3));
+        xGradient = sqrt(sum(xGradient .* xGradient, 3));
+    end
+
+    % Compute polar-coordinate representation
+    gradientMagnitude = sqrt(yGradient .* yGradient + xGradient .* xGradient);
+    gradientOrientation = atan2(xGradient, yGradient);
+end