diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..ed6eb3d
--- /dev/null
+++ b/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.
\ No newline at end of file
diff --git a/SemanticSoftSegmentation.m b/SemanticSoftSegmentation.m
new file mode 100644
index 0000000..ed306b4
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/SpectralMatting.m b/SpectralMatting.m
new file mode 100644
index 0000000..1069380
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/Superpixels.m b/Superpixels.m
new file mode 100644
index 0000000..ba661b7
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/affinityMatrixToLaplacian.m b/affinityMatrixToLaplacian.m
new file mode 100644
index 0000000..8d2a64f
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/demo.m b/demo.m
new file mode 100644
index 0000000..9a8a774
--- /dev/null
+++ b/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');
diff --git a/docia.png b/docia.png
new file mode 100644
index 0000000..65beb76
Binary files /dev/null and b/docia.png differ
diff --git a/groupSegments.m b/groupSegments.m
new file mode 100644
index 0000000..518bb6a
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/imageGradient.m b/imageGradient.m
new file mode 100644
index 0000000..094d40a
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/mattingAffinity.m b/mattingAffinity.m
new file mode 100644
index 0000000..d08bdb4
--- /dev/null
+++ b/mattingAffinity.m
@@ -0,0 +1,101 @@
+
+% Matting Affinity
+% This function implements the image matting approach described in
+% Anat Levin, Dani Lischinski, Yair Weiss, "A Closed Form Solution to 
+% Natural Image Matting", IEEE TPAMI, 2008.
+% All parameters other than image are optional. The output is a sparse
+% matrix which has non-zero element for the non-local neighbors of
+% the pixels given by binary map inMap.
+% - windowRadius defines the size of the window where the local normal 
+%   distributions are estimated.
+% - epsilon defines the regularization coefficient used before inverting
+%   covariance matrices. It should be larger for noisy images.
+
+function W = mattingAffinity(image, inMap, windowRadius, epsilon)
+
+    if ~exist('windowRadius', 'var') || isempty(windowRadius)
+        windowRadius = 1;
+    end
+    if ~exist('epsilon', 'var') || isempty(epsilon)
+        epsilon = 1e-7;
+    end
+    image = im2double(image);
+    windowSize = 2 * windowRadius + 1;
+    neighSize = windowSize^2;
+    [h, w, c] = size(image);
+    N = h * w;
+    epsilon = epsilon / neighSize;
+    
+    % No need to compute affinities in known regions if a trimap is defined
+    if nargin < 2 || isempty(inMap)
+        inMap = true(size(image, 1), size(image, 2));
+    end
+    
+    [meanImage, covarMat] = localRGBnormalDistributions(image, windowRadius, epsilon);
+
+    % Determine pixels and their local neighbors
+    indices = reshape((1 : h * w), [h w]);
+    neighInd = im2col(indices, [windowSize windowSize], 'sliding')';
+    inMap = inMap(windowRadius + 1 : end - windowRadius, windowRadius + 1 : end - windowRadius);
+    neighInd = neighInd(inMap, :);
+    inInd = neighInd(:, (neighSize + 1) / 2);
+    pixCnt = size(inInd, 1);
+
+    % Prepare in & out data
+    image = reshape(image, [N, c]);
+    meanImage = reshape(meanImage, [N, c]);
+    flowRows = zeros(neighSize, neighSize, pixCnt);
+    flowCols = zeros(neighSize, neighSize, pixCnt);
+    flows = zeros(neighSize, neighSize, pixCnt);
+
+    % Compute matting affinity
+    for i = 1 : size(inInd, 1)
+        neighs = neighInd(i, :);
+        shiftedWinColors = image(neighs, :) - repmat(meanImage(inInd(i), :), [size(neighs, 2), 1]);
+        flows(:, :, i) = shiftedWinColors * (covarMat(:, :, inInd(i)) \ shiftedWinColors');
+        neighs = repmat(neighs, [size(neighs, 2), 1]);
+        flowRows(:, :, i) = neighs;
+        flowCols(:, :, i) = neighs';
+    end
+    flows = (flows + 1) / neighSize;
+    
+    W = sparse(flowRows(:), flowCols(:), flows(:), N, N);
+    % Make sure it's symmetric
+    W = (W + W') / 2;
+end
+
+
+function [meanImage, covarMat] = localRGBnormalDistributions(image, windowRadius, epsilon)
+
+    if ~exist('windowRadius', 'var') || isempty(windowRadius)
+        windowRadius = 1;
+    end
+    if ~exist('epsilon', 'var') || isempty(epsilon)
+        epsilon = 1e-8;
+    end
+
+    [h, w, ~] = size(image);
+    N = h * w;
+    windowSize = 2 * windowRadius + 1;
+
+    meanImage = imboxfilt(image, windowSize);
+    covarMat = zeros(3, 3, N);
+
+    for r = 1 : 3
+        for c = r : 3
+            temp = imboxfilt(image(:, :, r).*image(:, :, c), windowSize) - meanImage(:,:,r) .*  meanImage(:,:,c);
+            covarMat(r, c, :) = temp(:);
+        end
+    end
+
+    for i = 1 : 3
+        covarMat(i, i, :) = covarMat(i, i, :) + epsilon;
+    end
+
+    for r = 2 : 3
+        for c = 1 : r - 1
+            covarMat(r, c, :) = covarMat(c, r, :);
+        end
+    end
+
+end
\ No newline at end of file
diff --git a/preprocessFeatures.m b/preprocessFeatures.m
new file mode 100644
index 0000000..09d670e
--- /dev/null
+++ b/preprocessFeatures.m
@@ -0,0 +1,40 @@
+
+% Pre-processing hyper-dimensional semantic feature vectors
+% as described in Section 3.5.
+
+function simp = preprocessFeatures(features, image)
+    % Filter out super high numbers due to some instability in the network
+    features(features < -5) = -5;
+    features(features > 5) = 5;
+    
+    % Filter each channel of features with image as the guide
+    if exist('image', 'var') && ~isempty(image)
+        fd = size(features, 3);
+        maxfd = fd - rem(fd, 3);
+        for i = 1 : 3 : maxfd
+            features(:, :, i : i+2) = imguidedfilter(features(:, :, i : i+2), image, 'NeighborhoodSize', 10);
+        end
+        for i = maxfd + 1 : fd
+            features(:, :, i) = imguidedfilter(features(:, :, i), image, 'NeighborhoodSize', 10);
+        end
+    end
+
+    % Run PCA and normalize to [0, 1]
+    simp = featuresPCA(features, 3);
+    for i = 1 : 3
+        simp(:,:,i) = simp(:,:,i) - min(min(simp(:,:,i)));
+        simp(:,:,i) = simp(:,:,i) / max(max(simp(:,:,i)));
+    end
+end
+
+function pcafeat = featuresPCA(features, dim)
+    features = double(features);
+    [h, w, d] = size(features);
+    features = reshape(features, [h*w, d]);
+    featmean = mean(features, 1);
+    features = features - ones(h*w, 1) * featmean;
+    covar = features' * features;
+    [eigvecs, ~] = eigs(covar, dim, 'LA');
+    pcafeat = features * eigvecs;
+    pcafeat = reshape(pcafeat, [h, w, dim]);
+end
diff --git a/readme.md b/readme.md
new file mode 100644
index 0000000..8979e1c
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,36 @@
+Semantic Soft Segmentation
+=================================================================
+
+This repository includes the spectral segmentation approach presented in
+
+    Yagiz Aksoy, Tae-Hyun Oh, Sylvain Paris, Marc Pollefeys and Wojciech Matusik, "Semantic Soft Segmentation", ACM Transactions on Graphics (Proc. SIGGRAPH), 2018 
+
+The network for semantic feature generation can be found [[here](https://github.com/iyah4888/SIGGRAPH18SSS)].
+
+Please refer to the [[project page](http://people.inf.ethz.ch/aksoyy/sss/)] for more information and example data.
+
+The Superpixels class requires [[ImageGraphs](http://www.mathworks.com/matlabcentral/fileexchange/53614-image-graphs)].
+
+License and citation
+------------
+
+This toolbox is provided for academic use only.
+If you use this code, please cite our paper:
+
+    @ARTICLE{sss,
+    author={Ya\u{g}{\i}z Aksoy and Tae-Hyun Oh and Sylvain Paris and Marc Pollefeys and Wojciech Matusik},
+    title={Semantic Soft Segmentation},
+    journal={ACM Transactions on Graphics (Proc. SIGGRAPH)},
+    year={2018},
+    pages = {72:1-72:13},
+    volume = {37},
+    number = {4}
+    }
+
+Credit
+------------
+Parts of this implementation are taken from
+
+    Anat Levin, Alex Rav-Acha, Dani Lischinski, "Spectral Matting", IEEE TPAMI, 2008
+
+The original source code for Spectral Matting can be found [[here](http://www.vision.huji.ac.il/SpectralMatting/)].
\ No newline at end of file
diff --git a/softSegmentsFromEigs.m b/softSegmentsFromEigs.m
new file mode 100644
index 0000000..92fcef4
--- /dev/null
+++ b/softSegmentsFromEigs.m
@@ -0,0 +1,187 @@
+
+% This function implements the constrained sparsification descibed in Section 3.4.
+% This methodology is intriduced by Levin et al. in "Spectral Matting", and
+% this function is an adapted version of their original source code. Find it here:
+% http://www.vision.huji.ac.il/SpectralMatting/
+
+function [softSegments, initialSegments] = softSegmentsFromEigs(eigVecs, eigVals, Laplacian, h, w, features, compCnt, maxIter, sparsityParam, imageGrad, initialSegments)
+    if (~exist('maxIter','var') || isempty(maxIter))
+      maxIter = 20;
+    end
+    if (~exist('sparsityParam','var') || isempty(sparsityParam))
+        sparsityParam = 0.9;
+    end
+    if (~exist('compCnt','var') || isempty(compCnt))
+        compCnt = 10;
+    end
+    eigValCnt = size(eigVecs, 2);
+    if (~exist('eigVals','var') || isempty(eigVals))
+        eigVals = 1e-5*eye(eigValCnt);
+    end
+
+    % Initialize with k-means if an initialization (initialSegments) isn't provided
+    if ~exist('initialSegments', 'var') || isempty(initialSegments)
+      % Use the features for the k-means, if provided
+      if ~(~exist('features', 'var') || isempty(features))
+        disp('          Computing k-means initialization using semantic features...');
+        %features = featuresPCA(features, 10);
+        features = reshape(features, [size(features, 1) * size(features, 2), size(features, 3)]);
+        [initialSegments, ~, ~, ~] = fastKmeans(features, compCnt);
+      else
+      % This is the default initialization from Spectral Matting
+        disp('          Computing k-means initialization...');
+        eigVals = abs(eigVals); % Sometimes there are -epsilon numbers that mess up the thing
+        initEigsCnt = 20;
+        initEigsWeights = diag(1 ./ diag(eigVals(2 : initEigsCnt + 1, 2 : initEigsCnt + 1) .^ 0.5 ));
+        initEigVecs = eigVecs(:, 2 : initEigsCnt + 1) * initEigsWeights;
+        [initialSegments, ~, ~, ~] = fastKmeans(initEigVecs, compCnt);
+      end
+    end
+    softSegments = zeros(length(initialSegments), compCnt);
+    for i = 1 : compCnt
+        softSegments(:, i) = double(initialSegments == i);
+    end
+
+    if ~exist('imageGrad', 'var') || isempty(imageGrad)
+        spMat = sparsityParam;
+    else
+        imageGrad(imageGrad > 0.2) = 0.2;
+        imageGrad = imageGrad + 0.8;
+        spMat = repmat(imageGrad(:), [1, compCnt]);
+    end
+
+    % Second derivative of first and second terms in the sparsity penalty
+    disp('          Starting optimization...');
+    thr_e = 1e-10;
+    w1 = 0.3;
+    w0 = 0.3;
+    e1 = w1 .^ sparsityParam * max(abs(softSegments-1), thr_e) .^ (spMat - 2);
+    e0 = w0 .^ sparsityParam * max(abs(softSegments), thr_e) .^ (spMat - 2);
+
+    scld = 1;
+    eig_vectors = eigVecs(:, 1 : eigValCnt);
+    eig_values = eigVals(1 : eigValCnt, 1 : eigValCnt);
+
+    % First iter no for removing zero components
+    removeIter = ceil(maxIter / 4);
+    removeIterCycle = ceil(maxIter / 4);
+
+    % Compute matting component with sparsity prior
+    for iter = 1 : maxIter
+        if rem(iter, 10) == 0
+            disp(['               Iteration ' int2str(iter) ' of ' int2str(maxIter)]);
+        end
+        % Construct the matrices in Eq(9) in Spectral Matting
+        tA = zeros((compCnt - 1) * eigValCnt);
+        tb = zeros((compCnt - 1) * eigValCnt, 1);
+        for k = 1 : compCnt - 1
+            weighted_eigs = repmat(e1(:, k) + e0(:, k), 1, eigValCnt) .* eig_vectors;
+            tA((k-1) * eigValCnt + 1 : k * eigValCnt, (k-1) * eigValCnt + 1 : k * eigValCnt) = eig_vectors' * weighted_eigs + scld * eig_values;
+            tb((k-1) * eigValCnt + 1 : k * eigValCnt) = eig_vectors' * e1(:,k);
+        end 
+        k = compCnt;
+        weighted_eigs = repmat(e1(:, k) + e0(:, k), 1, eigValCnt) .* eig_vectors;
+        ttA = eig_vectors' * weighted_eigs + scld * eig_values;
+        ttb = eig_vectors' * e0(:, k) + scld * sum(eig_vectors' * Laplacian, 2);
+
+        tA = tA + repmat(ttA, [compCnt - 1, compCnt - 1]);
+        tb = tb + repmat(ttb, [compCnt - 1, 1]);
+
+        % Solve for weights
+        y = reshape(tA \ tb, eigValCnt, compCnt - 1);
+
+        % Compute the matting comps from weights
+        softSegments = eigVecs(:, 1 : eigValCnt) * y;
+        softSegments(:, compCnt) = 1 - sum(softSegments(:, 1 : compCnt - 1), 2); % Sets the last one as 1-sum(others), guaranteeing \sum(all) = 1
+
+        % Remove matting components which are close to zero every once in a while
+       if iter > removeIter
+           nzii = find(max(abs(softSegments)) > 0.1);
+           compCnt = length(nzii);
+          softSegments = softSegments(:, nzii);
+           removeIter = removeIter + removeIterCycle;
+       end
+
+        % Recompute the derivatives of sparsity penalties
+        if length(spMat(:)) == 1
+          e1 = w1 .^ sparsityParam * max(abs(softSegments-1), thr_e) .^ (spMat - 2);
+          e0 = w0 .^ sparsityParam * max(abs(softSegments), thr_e) .^ (spMat - 2);
+        else
+          e1 = w1 .^ sparsityParam * max(abs(softSegments-1), thr_e) .^ (spMat(:, 1 : size(softSegments, 2)) - 2);
+          e0 = w0 .^ sparsityParam * max(abs(softSegments), thr_e) .^ (spMat(:, 1 : size(softSegments, 2)) - 2);
+        end
+    end
+    softSegments = reshape(softSegments, [h, w, size(softSegments, 2)]);
+end
+
+function [idx, C, sumd, D]=fastKmeans(X,K)
+  % X: points in the N-by-P data matrix
+  % idx - an N-by-1 vector containing the cluster indices of each point
+  %
+  % c - the K cluster centroid locations in a K-by-P matrix.
+  % sumd - the within-cluster sums of point-to-centroid distances in a 1-by-K vector.
+  % distances from each point to every centroid in a N-by-K.
+
+  startK = 5;
+  startK = min(K,startK);
+  maxIters = 100;  % Defualt of matlab is 100    
+
+    X = sign(real(X)) .* abs(X);
+  
+  [idx, C, sumd, D]=kmeans(X,startK,'EmptyAction','singleton','Start','cluster', ...
+      'Maxiter', maxIters, 'Replicates', 7);
+
+  valid_vec = zeros(1,startK);
+  scr_vec = zeros(1,startK)-1;
+
+  for compCnt = startK+1:K
+      % create a new cluster by splitting each cluster to two...
+      max_scr=-1;
+      clear min_C;
+      for cl = 1:compCnt-1
+        cl_mask = idx == cl;
+        cl_idxs = find(cl_mask);
+        clX = X(cl_idxs,:);
+        if (size(clX,1)> 2*size(clX,2))
+          if (valid_vec(cl) == 0)
+            [tmp_idx, tmp_C, ~, tmp_D]=kmeans(clX,2,'EmptyAction','singleton','Start','cluster', 'Maxiter', maxIters);
+            % chk how much the partition helps ...
+            scr=sum(min(D(cl_idxs,:),[],2))-sum(min(tmp_D,[],2));
+            scr_vec(cl) = scr;
+          else % we already saved it...
+            scr = scr_vec(cl);
+          end         
+        else
+          scr=-2;
+          scr_vec(cl) = scr;
+        end  
+        if (scr > max_scr)
+          if (valid_vec(cl)==1) % not for the scr. Just for the idxs.
+            [tmp_idx, tmp_C, ~, ~]=kmeans(clX,2,'EmptyAction','singleton','Start','cluster', 'Maxiter', maxIters);
+          end
+          
+          max_scr = scr;
+          bestC = [C;tmp_C(2,:)];  
+          bestC(cl,:) = tmp_C(1,:);
+          best_cl = cl;         
+
+          best_idx = idx;
+          best_idx(cl_idxs) = (tmp_idx == 1)*best_cl + (tmp_idx == 2)*compCnt;
+        end
+        valid_vec(cl) = 1;
+      end        
+      C = bestC;
+      idx = best_idx;
+      
+      valid_vec = [valid_vec, 0];   % the two new clusers are new, so their
+      valid_vec(best_cl) = 0;      % score have not been computed yet.
+      scr_vec = [scr_vec, -1];
+      scr_vec(best_cl) = -1;
+
+      if (compCnt < 13)
+        [idx, C, sumd, D]=kmeans(X, compCnt, 'EmptyAction', 'singleton', 'Start', C, 'Maxiter', maxIters);       
+        valid_vec = zeros(1,compCnt);
+      end
+  end
+
+end
\ No newline at end of file
diff --git a/sparsifySegments.m b/sparsifySegments.m
new file mode 100644
index 0000000..92844b2
--- /dev/null
+++ b/sparsifySegments.m
@@ -0,0 +1,69 @@
+
+% This function implements the relaxed sparsification descibed in Section 3.4
+
+function softSegments = sparsifySegments(softSegments, Laplacian, imageGrad)
+ 
+    sigmaS = 1; % sparsity
+    sigmaF = 1; % fidelity
+    delta = 100; % constraint
+    [h, w, compCnt] = size(softSegments);
+    N = h * w * compCnt;
+
+    if ~exist('imageGrad', 'var') || isempty(imageGrad)
+        % If image gradient is not provided, set the param to the default 0.9
+        spPow = 0.90;
+    else
+        % Compute the per-pixel sparsity parameter from the gradient
+        imageGrad(imageGrad > 0.1) = 0.1;
+        imageGrad = imageGrad + 0.9;
+        spPow = repmat(imageGrad(:), [compCnt, 1]);
+    end
+
+    % Iter count for pcg and main optimization
+    itersBetweenUpdate = 100;
+    highLevelIters = 20;
+
+    % Get rid of very low/high alpha values and normalize
+    softSegments(softSegments < 0.1) = 0;
+    softSegments(softSegments > 0.9) = 1;
+    softSegments = softSegments ./ repmat(sum(softSegments, 3), [1 1 size(softSegments, 3)]);
+    
+    % Construct the linear system
+    lap = Laplacian;
+    for i = 2 : compCnt
+        Laplacian = blkdiag(Laplacian, lap);
+    end
+
+    % The alpha constraint
+    C = repmat(speye(h*w), [1 compCnt]);
+    C = C' * C;
+    Laplacian = Laplacian + delta * C;
+
+    % The sparsification optimization
+    softSegments = softSegments(:);
+    compInit = softSegments; % needed for fidelity energy
+    for iter = 1 : highLevelIters
+        if rem(iter, 5) == 0
+            disp(['               Iteration ' int2str(iter) ' of ' int2str(highLevelIters)]);
+        end
+        [u, v] = getUandV(softSegments, spPow); % The sparsity energy
+        A = Laplacian + sigmaS * (spdiags(u, 0, N, N) + spdiags(v, 0, N, N)) + sigmaF * speye(N);
+        b = sigmaS * v + sigmaF * compInit + delta;
+        [softSegments, ~] = pcg(A, b, [], itersBetweenUpdate, [], [], softSegments);
+    end
+
+    % One final iter for good times (everything but sparsity)
+    A = Laplacian + sigmaF * speye(N);
+    b = sigmaF * softSegments + delta;
+    softSegments = pcg(A, b, [], 10 * itersBetweenUpdate, [], [], softSegments);
+
+    % Ta-dah
+    softSegments = reshape(softSegments, [h w compCnt]);
+end
+
+function [u, v] = getUandV(comp, spPow)
+    % Sparsity terms in the energy
+    eps = 1e-2;
+    u = max(abs(comp(:)), eps) .^ (spPow - 2);
+    v = max(abs(1 - comp(:)), eps) .^ (spPow - 2);
+end
\ No newline at end of file
diff --git a/visualizeEigenvectorRedGreen.m b/visualizeEigenvectorRedGreen.m
new file mode 100644
index 0000000..718122e
--- /dev/null
+++ b/visualizeEigenvectorRedGreen.m
@@ -0,0 +1,9 @@
+
+function vis = visualizeEigenvectorRedGreen(eigenvector)
+    % negative values are shown in red, and positive in green
+    vis = zeros(size(eigenvector, 1), size(eigenvector, 2), 3);
+    vis(:,:,1) = -100 * eigenvector;
+    vis(:,:,2) = 100 * eigenvector;
+    vis(vis < 0) = 0;
+    vis(vis > 1) = 1;
+end
\ No newline at end of file
diff --git a/visualizeSoftSegments.m b/visualizeSoftSegments.m
new file mode 100644
index 0000000..5f14e84
--- /dev/null
+++ b/visualizeSoftSegments.m
@@ -0,0 +1,43 @@
+
+% Assigns a distinct solid color to eact soft segment and composites
+% them using the corresponding alpha calues
+
+function [vis, softSegments] = visualizeSoftSegments(softSegments, doOrdering)
+    if ~exist('doOrdering', 'var') || isempty(doOrdering)
+        doOrdering = false;
+    end
+    if doOrdering
+        % Order layers w.r.t. sum(alpha(:)) -- makes visualizations more consistent across images
+        order = determineAlphaOrder(softSegments);
+        softSegments = softSegments(:, :, order);
+    end
+    % A hard-coded set of 'distinct' colors
+    colors = createPalette();
+    % One solid color per segment, mixed w.r.t. alpha values 
+    vis = repmat(softSegments(:,:,1), [1 1 3]) .* repmat(colors(1,1,:), [size(softSegments, 1), size(softSegments, 2), 1]);
+    for i = 2 : size(softSegments, 3)
+        vis = vis + repmat(softSegments(:,:,i), [1 1 3]) .* repmat(colors(i,1,:), [size(softSegments, 1), size(softSegments, 2), 1]);    
+    end
+end
+
+function colors = createPalette()
+    % https://graphicdesign.stackexchange.com/questions/3682/where-can-i-find-a-large-palette-set-of-contrasting-colors-for-coloring-many-d
+    % From P. Green-Armytage (2010): A Colour Alphabet and the Limits of Colour Coding. 
+    % Colour: Design & Creativity (5) (2010): 10, 1-23
+    % http://eleanormaclure.files.wordpress.com/2011/03/colour-coding.pdf
+    colors = reshape(... 
+                                    [0,117,220; 255,255,128; 43,206,72; 153,0,0; 128,128,128; 240,163,255; 153,63,0; 76,0,92;...
+                                    0,92,49; 255,204,153; 148,255,181;...
+                                    143,124,0; 157,204,0; 194,0,136; 0,51,128; 255,164,5;...
+                                    255,168,187; 66,102,0; 255,0,16; 94,241,242; 0,153,143;...
+                                    224,255,102; 116,10,255; 255,255,0; 255,80,5; 25,25,25]...
+                                    , [26, 1, 3]) / 255;
+    colors = [colors; colors];
+end
+
+function order = determineAlphaOrder(alphas)
+    alphas = sum(alphas, 1);
+    alphas = sum(alphas, 2);
+    alphas = reshape(alphas, [size(alphas, 3) 1]);
+    [~, order] = sort(alphas, 'descend');
+end
\ No newline at end of file