From f618466c25d43f3bae9e40920273bf77de1e1149 Mon Sep 17 00:00:00 2001
From: leochanj105 <leochanj@live.unc.edu>
Date: Mon, 19 Oct 2020 23:09:30 -0400
Subject: initial sd-vbs

initial sd-vbs

add sd-vbs

sd-vbs
---
 SD-VBS/common/toolbox/MultiNcut/computeMultiW.m | 245 ++++++++++++++++++++++++
 1 file changed, 245 insertions(+)
 create mode 100755 SD-VBS/common/toolbox/MultiNcut/computeMultiW.m

(limited to 'SD-VBS/common/toolbox/MultiNcut/computeMultiW.m')

diff --git a/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m b/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m
new file mode 100755
index 0000000..b2cb7ea
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/computeMultiW.m
@@ -0,0 +1,245 @@
+%inputs : image, number of layers, distance defining the subgrid, the edge filter scales for each layer, radius for each layer,
+%edge variance for filter, shape of the neighbourhood layout ('square', 'star', 'cross'), sigma for intensity affinity,
+% sigma for distance influence in affinity, weight coefficients for Wpps in the multiscale matrix.
+%output : multiscale affinity matrix , extern constraint matrix, affinity matrices of each layer seperatly.
+
+
+function [multiWpp,constraintMat, Wind,data,emag,ephase]= computeMultiW (image,data);
+
+%variables
+
+if  isempty(data.layers.number)
+    n=2;
+else
+    n=data.layers.number;
+end
+
+if isempty(data.layers.dist)
+    dist_grid=3;
+else
+    dist_grid=data.layers.dist;
+end
+
+if isempty(data.W.scales)
+    s=1:n;
+elseif (length(data.W.scales)==n)
+    s=data.W.scales;
+else
+    s=1:n;
+end
+
+if isempty(data.W.radius)
+    r(1)=2;
+    for i=2:n
+        r(i)=10;
+    end
+else
+    r=data.W.radius;
+end
+
+
+if isempty(data.W.edgeVariance)
+    data.W.edgeVariance=0.1;
+end
+
+if isempty(data.W.gridtype)
+    data.W.gridtype='square';
+end
+
+if isempty(data.W.sigmaI)
+    data.W.sigmaI=0.12;
+end
+
+if isempty(data.W.sigmaX)
+    data.W.sigmaX=10;
+end
+
+if isempty(data.layers.weight)
+    coef(1)=5;
+    coef(2:n)=200;
+elseif (length(data.layers.weight)==n)
+    coef=data.layers.weight;
+else
+    coef(1)=5;
+    coef(2:n)=100;  %200  
+end
+
+if isempty(data.W.mode)
+    data.W.mode=mixed;
+end
+
+
+[p1,q1,ignore]=size(image);
+image=image(:,:,1);
+filter_par = [4,30,4];  %[9,30,4]
+[x,y,gx,gy,par,threshold,emag,ephase,g,FIe,FIo,mago] = quadedgep2(image,filter_par,data,0.001);
+minW=10^(-2); %-3
+
+
+% function [multiWpp,constraintMat,p,q,Wppp,subgrid] = computemultiWpp(n,imageX,r,dist_grid,s,dataWpp,emag,ephase,minW,mode,facteurMul,contrainte,tt,gridtype,colormode,imageOriginale,subgridImageReduite,pG,qG)
+
+p= p1*ones(n,1);
+q= q1*ones(n,1);
+d= dist_grid*ones(n,1);
+d(1)=1;
+for (i=2:n)
+    d(i)=d(i)*3^(i-2);
+end
+p=ceil(p1./d);
+q=ceil(q1./d);
+
+%computation of the subgrids (the first pixel is coded by one). S{i,j}(k) is the index of
+%the kth pixel of the jth grid in the ith grid.
+
+for i=1:n-1
+    for j=i+1:n
+        a=[0:p(j)*q(j)-1];
+        subgrid{i,j}=p(i)*(floor(a/p(j)))*d(j)/d(i)+(1+mod(a,p(j))*d(j)/d(i));
+    end
+end
+
+%computation of the independent W matrix for each layer Wind{i} 1=<i=<n.
+
+[w1i,w1j]=cimgnbmap([p1,q1], r(1), 1);
+
+if strcmp(data.W.mode,'mixed')
+    rMin = 0;
+    imageXX=double(image(:));
+    sigmaI= (std(imageXX(:)) + 1e-10 )* data.W.sigmaI;
+    Wpp{1} = multiIntensityFirstLayer(double(image),w1i,w1j,rMin,data.W.sigmaX,sigmaI,minW);
+    Wpp2= affinityic(emag(:,:,s(1)),ephase(:,:,s(1)),w1i,w1j,max(max(emag(:,:,s(1)))) * data.W.edgeVariance);
+    Wpp{1} = sqrt(Wpp{1} .* Wpp2)+0.1*Wpp2;
+    
+elseif strcmp(data.W.mode,'notmixed')
+    Wpp{1}= affinityic(emag(:,:,s(1)),ephase(:,:,s(1)),w1i,w1j,max(max(emag(:,:,s(1)))) * data.W.edgeVariance);
+    
+elseif strcmp(data.W.mode,'intensity')
+    rMin = 0;
+    imageXX=double(image(:));
+    sigmaI= (std(imageXX(:)) + 1e-10 )* data.W.sigmaI;
+    Wpp{1} = multiIntensityFirstLayer(double(image),w1i,w1j,rMin,data.W.sigmaX,sigmaI,minW);  
+    
+end
+Wpp{1}=coef(1)*(Wpp{1}+Wpp{1}')/2;
+%Wpp{1}= coef(1)*Wpp{1};
+Wind{1}=Wpp{1};
+
+
+
+
+for i=2:n 
+    if strcmp(data.W.gridtype,'square')
+        [wwi,wwj]=cimgnbmap([p(i),q(i)], r(i), 1); 
+    elseif strcmp(data.W.gridtype,'star')
+        [wwi,wwj]=cimgnbmap_star([p(i),q(i)], r(i), 1);
+    elseif strcmp(data.W.gridtype,'cross')
+        [wwi,wwj]=cimgnbmap_cross([p(i),q(i)], r(i), 1);
+    end
+    wwi=double(wwi);
+    wiInOriginalImage=(p1*(floor(wwi/p(i)))*d(i))+(1+mod(wwi,p(i))*d(i));
+        wiInOriginalImage=(p1*(floor(wwi/p(i)))*d(i))+(1+mod(wwi,p(i))*d(i));
+
+    wiInOriginalImage= uint32(wiInOriginalImage);
+
+if strcmp(data.W.mode,'mixed')
+    Wpp2= multiAffinityic(emag(:,:,i),ephase(:,:,i),wiInOriginalImage,wwj,subgrid{1,i},p(i),q(i),uint32(wwi),max(max(emag(:,:,i))) * data.W.edgeVariance);
+    a=floor(d(i)/d(i-1));
+    if (mod(a,2)==0)
+        a=a+1;
+    end
+%     Wpp{i} = multiIntensityWppc(double(imageX),wiInOriginalImage,wwj,rMin,dataWpp.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),wi{i});
+    
+    Wpp{i} = multiIntensityWppc(double(image),wiInOriginalImage,wwj,rMin,data.W.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),uint32(wwi));
+    Wpp{i} = sqrt(Wpp{i} .* Wpp2)+0.1*Wpp2;
+elseif strcmp(data.W.mode,'notmixed')
+    Wpp{i}= multiAffinityic(emag(:,:,i),ephase(:,:,i),wiInOriginalImage,wwj,subgrid{1,i},p(i),q(i),uint32(wwi),max(max(emag(:,:,i))) * data.W.edgeVariance);
+elseif strcmp(data.W.mode,'intensity')
+    Wpp{i} = multiIntensityWppc(double(image),wiInOriginalImage,wwj,rMin,data.W.sigmaX,sigmaI,minW,subgrid{1,i},p(i),q(i),uint32(wwi));
+end
+Wpp{i}= coef(i)*(Wpp{i}+Wpp{i}')/2;
+
+%Wpp{i}= coef(i)*Wpp{i};
+Wind{i}=Wpp{i};
+
+end
+
+%computation of the intern contraint matrices C{i,j}.
+
+for i=1:n-1
+    r=floor(d(i+1)/(d(i)*2));
+    [wwi,wwj]=cimgnbmap([p(i),q(i)], r, 1);
+    wi{i}=wwi;
+    wj{i}=wwj;
+end
+
+for i=1:n-1
+    for j=i+1:n
+        C{i,j}=sparse(p(i)*q(i),p(j)*q(j));
+        firstPointer=double(wj{i}(subgrid{i,j}))+1;   
+        lastPointer=double(wj{i}(subgrid{i,j}+1));
+        invNbNeighbours=1./(lastPointer-firstPointer+1);
+        for (k=1:p(j)*q(j))
+            for (m=firstPointer(k):lastPointer(k))
+                C{i,j}(double(wi{i}(m))+1,k)=invNbNeighbours(k);
+            end
+        end
+    end
+end
+            
+%Assembling the built matrices to make up multiWpp.
+for i=1:n
+    if (i>1)
+        for j=i-1:-1:1
+            Wpp{i}=[C{j,i}',Wpp{i}];
+        end
+    end
+    if (i<n)
+        for j=i+1:n
+            Wpp{i}=[Wpp{i},C{i,j}];
+        end
+    end
+end
+
+% %Assembling the built matrices to make up Wpp without intern constrains.
+% for i=1:n
+%     if (i>1)
+%         for j=i-1:-1:1
+%             Wpp{i}=[sparse(p(i)*q(i),p(j)*q(j)),Wpp{i}];
+%         end
+%     end
+%     if (i<n)
+%         for j=i+1:n
+%             Wpp{i}=[Wpp{i},sparse(p(i)*q(i),p(j)*q(j))];
+%         end
+%     end
+% end
+
+clear Wind;Wind = 1;
+
+multiWpp=Wpp{1}; clear Wpp{1}
+for i=2:n
+    multiWpp=[multiWpp;Wpp{i}];clear Wpp{i}
+end
+
+ 
+% Computing the average extern constraint
+
+ pq=sum(p(2:n).*q(2:n));
+ p2q2=p(2)*q(2);
+ constraintMat=[-C{1,2};speye(p2q2);sparse(pq-p2q2,p2q2)];
+ if n>2
+     for i=3:n
+         piqi=p(i)*q(i);
+         if i~=n
+             constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi);sparse(pq-sum(p(2:i).*q(2:i)),piqi)]];
+         else
+             constraintMat=[constraintMat,[sparse(sum(p(1:i-2).*q(1:i-2)),piqi);-C{i-1,i};speye(piqi)]];
+         end
+     end
+ end
+
+ % saving useful information
+ %subgrids, p and q
+ data.subgrid=subgrid;
+ data.p=p;
+ data.q=q;
-- 
cgit v1.2.2