expected_genmodel_fhedgecolsw.m

基于matlab写的对人体的分割程序

function [expWWW,expWCAll,expOrient,NN] = expected_genmodel_FHedgecolsW(im,genmodel)
%[expWWW,colIm,partIms,logNN] = expected_genmodel_FHedgecolsW(im,genmodel)
%Takes the RAW image, since it needs to compute both edges and olor
%Does both edge mask and a color mask
%Assumes genmodel has a allfgP and colmask field
%ALSO assumes genmodel has a .www feild, which is a scalar that weights the
%.ww masks; this is the only thing that is tweaked!


[imy imx imz] = size(im);

numClass = max(genmodel.equiv_class);
numTypes = size(genmodel.dag,1);
imo = size(genmodel.orient,2);
respIm = zeros([imy imx imo numTypes]);
assert(numClass == numTypes);
lens = genmodel.len*2; wids = genmodel.wid*2;

%Convert orient to exp scale (and normalize)
genmodel.orient = exp(genmodel.orient);
genmodel.orient = genmodel.orient./repmat(sum(genmodel.orient,2),1,imo);
%genmodel.orient(1,:) = 0; genmodel.orient(1,1) = 1;
%genmodel.orient(2,:) = 0; genmodel.orient(2,19) = 1;
%Flip the non-root ones since we'll be going convolution
%genmodel.orient = fliplr(genmodel.orient);
%genmodel.orient(genmodel.root,:) = fliplr(genmodel.orient(genmodel.root,:));
%genmodel.orient(2,:) = 1;

%Get out edge and color histogram map
%[m,t] = mydetGMc(im,2); m = double(imdilate(m > 10,ones(3)));
%imHist = imvq16(im);
%[dummy,dummy,imHist] = unique(imHist);
DEBUG = 0;
if DEBUG,
    [m,t] = mydetGMc(im,2); m = double(imdilate(m > 10/255,ones(3)));
    imHist = im+1;
else
    [m,t] = mydetGMc(im,2); m = double(imdilate(m > 10,ones(3)));
    imHist = imvq16(im);
    [dummy,dummy,imHist] = unique(imHist);
end

LOG_MAX = 500;
LOG_MIN = -500;

pNN = zeros(numTypes,1);
origWW = zeros(size(respIm));
for p = 1:numTypes,
    ww = genmodel.ww(:,:,p);
    mask = imresize(ww,2*[genmodel.len(p) 2*genmodel.wid(p)],'nearest');
    if any(mask(:)),
        mask = mask*(sum(abs(ww(:)))/sum(abs(mask(:))));
    end
    resp = getSegmentsEdge(m,mask);
    %Treat this as a feature that is weighted by genmodel.www
    origWW(:,:,:,p) = resp;
    resp = resp*genmodel.www(p);
    
    %Now take the color
    ww = genmodel.wc(:,:,p);
    mask = imresize(ww,2*[genmodel.len(p) 2*genmodel.wid(p)],'nearest');
    if any(mask(:)),
        mask = mask*(sum(abs(ww(:)))/sum(abs(mask(:))));
    end
    fgP = reshape(genmodel.fgP(imHist,find(genmodel.colmask(:,p))),imy,imx);
    resp = resp + getSegmentsEdge(fgP,mask);    
    resp = partshiftZ0(resp,genmodel.len(p),0); %Shift so rectangles are anchored at the top
    resp(resp > LOG_MAX) = LOG_MAX; resp(resp < LOG_MIN) = LOG_MIN;
    resp = exp(resp); 
    NN = sum(resp(:));
    pNN(p) = log(NN);
    respIm(:,:,:,p) = resp/NN;
end
%keyboard;
%tmp = respIm(:,:,:,2); [dummy,I] = max(tmp(:)); [y,x,a] = ind2sub([imy imx imo],I);
%respIm(:,:,:,p) = 0; respIm(100,40,24,p) = 1;
%respImOrig = respIm;
%imagesc(sum(respIm(:,:,:,p),3))

mean_y = round(genmodel.mu(:,2));
mean_x = round(genmodel.mu(:,1));

for p = genmodel.post(1:end-1),
  %Multiply in all the messages from the kids to the current respIm
  parType = genmodel.parents{p};
  %kids = find(genmodel.dag(p,:));
  %We sum over all orientations
  %currResp = sum(prod(respIm(:,:,:,[p kids]),4),3);
  %currResp = sum(respIm(:,:,:,p),3);
  currResp = reshape(respIm(:,:,:,p),imy*imx,imo);
  %Reflect borders
  %currResp = cat(2,currResp(:,imo/2:(end-1)),currResp,currResp(:,2:imo/2));
  currResp = cat(2,currResp(:,(imo/2+1):end),currResp,currResp(:,1:(imo/2-1)));

  %Convolve
  currResp = conv2(currResp,fliplr(genmodel.orient(p,:)),'valid');
  currResp = reshape(currResp,[imy imx imo]);  
    
  %Smooth it
  %std = sqrt(var(p));
  %kernal = exp(-1/(2*var(p))*[-3*std:3*std].^2); kernal = kernal/sum(kernal);
  %currResp = conv2(kernal',kernal,currResp,'same');
  currResp = local_sum_zero(currResp,genmodel.box(p)*2);
  
  %Shift it
  currResp = partshiftZ0(currResp,mean_y(p),mean_x(p));  
  %Normalize
  %currResp = currResp/(sum(currResp(:)) + eps);
  %Pass message to the parent
  respIm(:,:,:,parType) = respIm(:,:,:,parType).*currResp;
end

%Ignore root position and orientation
p = genmodel.root;
currResp = respIm(:,:,:,p);
%For now, zero out upside down roots
%currResp(:,:,[7:19]) = 0;
%Add in root potential
currResp = currResp.*repmat(permute(genmodel.orient(p,:),[1 3 2]),[imy imx 1]);

%currResp(:,:,[1:5 19:24]) = 0;
%currResp(:,:,[1:11 13:24]) = 0;%Debugging
%Save the normalizing constant
NN = sum(currResp(:));
respIm(:,:,:,p) = currResp/NN;
NN = log(NN) + sum(pNN);

%Pass messages down
expOrient = zeros(size(genmodel.orient));
expOrient(p,:) = sum(reshape(currResp,imy*imx,imo));
expOrient(p,:) = expOrient(p,:)/sum(expOrient(p,:));

%Define downstream orientations as walking counter-clockwise starting 
%from 13, pointing down
%ie [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]
%-> [1 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2]
%genmodel.orient = genmodel.orient(:,[1 end:-1:2]);
for p = genmodel.pre(2:end),
    
    parType = genmodel.parents{p};
    currResp = respIm(:,:,:,parType);
    currResp = partshiftZ0(currResp,-mean_y(p),-mean_x(p));
    currResp = reshape(currResp,imy*imx,imo);

    if nargout > 1,
        %Calculate the normalized child response
        kidResp = local_sum_zero(respIm(:,:,:,p),genmodel.box(p)*2);
        kidResp = reshape(kidResp,imy*imx,imo);
        %kidResp = cat(2,kidResp(:,imo/2:end-1),kidResp,kidResp(:,2:imo/2));
        kidResp = cat(2,kidResp(:,(imo/2+1):end),kidResp,kidResp(:,1:(imo/2-1)));
        normChild = conv2(kidResp,fliplr(genmodel.orient(p,:)),'valid') + eps;
        for i = 1:imo,
            aa = reshape(currResp,[imy imx imo]); bb =reshape(kidResp(:,i:(i+imo-1)),[imy imx imo]); cc = reshape(normChild,[imy imx imo]);
            tmp = currResp.*kidResp(:,i:(i+imo-1))./normChild;
            expOrient(p,i) = genmodel.orient(p,i)*sum(tmp(:));
        end
        %Flip back
        %expOrient(p,:) = expOrient(p,[1 end:-1:2]);
    end
    
    %Now do normal message passing
    %Circularly convolve
    %currResp = cat(2,currResp(:,imo/2:end-1),currResp,currResp(:,2:imo/2));
    currResp = cat(2,currResp(:,(imo/2+1):end),currResp,currResp(:,1:(imo/2-1)));
    %Shift so previous center (13) is now at center(1)
    currOrient = genmodel.orient(p,[1 end:-1:2]);
    currResp = conv2(currResp,fliplr(currOrient),'valid');
    %sum(currResp); %Useful debugging
    currResp = reshape(currResp,[imy imx imo]);
    currResp = local_sum_zero(currResp,genmodel.box(p)*2);
       
    %Pass it equally to all rotated copies
    currResp = respIm(:,:,:,p).*currResp;
    
    %Save the marginal
    respIm(:,:,:,p) = currResp./(sum(currResp(:))+eps);
end

%Shift the respIm to align with center of segment
for p = 1:numTypes,
    currResp = partshiftZ0(respIm(:,:,:,p),-genmodel.len(p),0);
    currResp = currResp./sum(currResp(:));
    respIm(:,:,:,p) = currResp;
end

if nargout == 1,
    partIms = zeros(imy,imx,numTypes);
    for p = 1:numTypes,
        len = genmodel.len(p);
        wid = genmodel.wid(p);
        %Condense into 180 degrees
        tmp = respIm(:,:,1:imo/2,p) + respIm(:,:,imo/2 + 1:end,p);
        kernal = ones(len*2,wid*2); %kernal = kernal/sum(kernal(:));
        %Spread by convolution
        for i = 1:imo/2,
            partIms(:,:,p) = partIms(:,:,p) + filter2(imrotate(kernal,(i-1)*360/imo),tmp(:,:,i));
        end
    end
    rgbs = ones(numTypes,1);
    ind = find(genmodel.dag(genmodel.root,:));
    rgbs(ind) = 2;
    ind = find(any(genmodel.dag(ind,:)));
    rgbs(ind) = 3;
    colIm = zeros(imy,imx);
    for i = 1:3,
        colIm(:,:,i) = sum(partIms(:,:,rgbs == i),3);
    end
    %expWWAll = colIm;
    %colIm = partIms;
    clear expWWW;
    expWWW.a = colIm;
    expWWW.b = partIms;
    expWWW.NN = NN;
    return;
end

expWWW = zeros(size(genmodel.www));
expWCAll = zeros(size(genmodel.wc));
for p = 1:numTypes,
    %p
    %expWWAll(:,:,p) = avgSegmentsEdge(m,lens(p),wids(p),respIm(:,:,:,p));
    tmp1 = respIm(:,:,:,p); tmp2 = origWW(:,:,:,p);
    expWWW(p) = tmp1(:)'*tmp2(:);
    fgP = reshape(genmodel.fgP(imHist,find(genmodel.colmask(:,p))),imy,imx);
    expWCAll(:,:,p) = avgSegmentsEdge(fgP,lens(p),wids(p),respIm(:,:,:,p));
end

DEBUG = 0;
if DEBUG,
    %Find the most likely pose
    [dummy,I] = max(currResp(:));
    [yy,xx,oo] = deal(zeros(numTypes,1));
    [yy(p),xx(p),oo(p)] = ind2sub(size(currResp),I);
    %keyboard;
    currRespOrig = currResp; 
    respImOrig = respIm;
    respIm = respImOrig;
    currResp = currRespOrig;
    for p = genmodel.pre(2:end),
        parType = genmodel.parents{p};
        currResp = respIm(:,:,:,p);
        %Contruct a kernal from the parent;
        kx = exp(-([1:imx] - xx(parType)).^2/(.5*var_x(p))); kx = kx/sum(kx);
        ky = exp(-([1:imy] - yy(parType)).^2/(.5*var_y(p))); ky = ky/sum(ky);

        
        paru = uu(oo(parType)); parv = vv(oo(parType));
        sq_do = 2*(1-(uu*paru + vv*parv)); 
        ko = exp(-sq_do./(.5*var_o(p))); ko = ko/sum(ko);
        kernal = repmat(ky(:),[1 length(kx) length(ko)]).*repmat(kx,[length(ky) 1 length(ko)]).*repmat(shiftdim(ko,-1),[length(ky) length(kx) 1]);
        currResp = currResp.*kernal;
        curr1 = currResp;
        %Translate back into image coordinates
        currResp = partshiftZ0(currResp,-mean_y(p),-mean_x(p));
        currShift = currResp;
        currResp = circshift(currResp,[0 0 mean_o(p)]);
        currCirc = currResp;
        [dummy,I] = max(currResp(:));
        [yy(p),xx(p),oo(p)] = ind2sub(size(currResp),I);
    end

    %Create segs
    seg.u = uu(oo)';
    seg.v = vv(oo)';
    seg.len = genmodel.len(genmodel.equiv_class);
    seg.w = genmodel.wid(genmodel.equiv_class);
    seg.x = xx + seg.u.*seg.len;
    seg.y = yy + seg.v.*seg.w;
    seg.type = [1:numTypes]';
    genmodel = seg;
end