2

I am very new with MRF and not that much good at programming. I have obtained probability map from semantic segmentation using a CNN, I have to optimize the segmentation by using Markov Random Fields (MRF). I download the code provided by Shai Bagon in this link GCmex. Energy minimization is performed based on either alpha expansion or swap.

I compiled the code by mex and I need to refine the Unary and pair-wise energy minimization functions. I have a stack of images and need to extract the 6-neighborhood grid and include the refined neighboring in the pair-wise function.

The input to the Unary function is the probability map which is a stack with size (256,256,4) for 4 different classes: enter image description here

My questions are: Has someone refined the code according to the defining different energy function 1) I wanna change Unary and pair-wise formulation). Which functions and which parts of code should be refined and recompiled again?

2) How to change the w_i,j? it is calculate based on intensity difference, here we have only probabilities, Is it the difference of probabilities of two adjacent voxels?

I really appreciate your help. Thanks

Shai
  • 111,146
  • 38
  • 238
  • 371
S.EB
  • 1,966
  • 4
  • 29
  • 54
  • How many 256x256x4 probability maps do you have? It seems like you have a map for each 2D slice, but you are doing 3D segmentation? how many slices do you have? Do you have the source images/slices for which you got the probability maps? – Shai Oct 24 '17 at 10:41
  • Hi Shai, thanks for responding, I am doing 2d segmentation by Cnn, which the output is 4 probability maps for one image (each probability map relates to one class), however, I have to do optimization by MRF using a grid of 6-neighbors (4 in current probability map and the other 2 in adjacent probability map of same class from the sequence. For the time being I want to do for each class separately (binary) to see how it works and the extend for all classes. How can I do it? which functions need to be changed for new data and smoothing terms?Thanks a lot. – S.EB Oct 24 '17 at 11:43
  • 6-connect is for 3D grids. Your grid is 2D, you can have either 4-connect or 8-connect. The graph only connects neighboring pixels and not neighboring entries in probability map. – Shai Oct 24 '17 at 11:47
  • It is a 6-neighbors in a grid of probability of a certain class (for example i). For example,If the current voxel is (r,c,k), the neighborhood includes 4 neighbors of the voxel plus two other neighbors from adjacent slices (probability of class i in slice k-1 and k+1, which are two voxel in (r,c,k+1) and (r,c,k-1). Since your GCmex code accepts an rgb image with discrete values, how can I refine it to take continuous probabilities [0,1] and which functions should change for refining potential terms? Thanks for your help – S.EB Oct 24 '17 at 11:53
  • so you have 4D probability map: `(r,c,k)x4`. do you have the input slices? what exactly do you have? – Shai Oct 24 '17 at 11:59
  • Yes, I have 60 images in a sequence and their ground truth (with 4 labels 0-3), each time I am sending a single image to CNN and output gives a stack of 4 probability maps per each class. Now I have to maximize the posterior prob. by using MRF (alpha-exp. or alpha-betha). All together, image, gt and stack of 4 as probability maps. – S.EB Oct 24 '17 at 12:12

1 Answers1

3

You have 60 slices of 256x256 pix (tot ~4G voxels), that is slices is a 256-by-256-by-60 array. Once you feed slices into your net (one by one or in batches - whatever works best for you) you have prob probability of size 256-by-256-by-60-by-4.
I suggest you use third constructor of GCMex to construct your graph for optimization.
To do so, you first need to define a sparse graph. Use sparse_adj_matrix:

[ii jj] = sparse_adj_matrix([256 256 60], 1, 1);  % 6-connect 3D grid
n = prod([256 256 60]);  % num voxels
wij = exp(-((slices(ii)-slices(jj)).^2)/(2*sig2));  % -|Ii-Ij|^2/2\sig^2
W = sparse(ii, jj, wij, n, n);  % sparse grid graph

Once you have the graph, it's all down hill from here:

Dc = -reallog(reshape(prob, n, 4)).';  %' unary/data term 
lambda = 2;  % relative weight of the smoothness term
gch = GraphCut('open', Dc, lambda*(ones(4)-eye(4)), W);  % construct the graph
[gch L] = GraphCut('expand', gch);  % minimize using "expand" method
gch = GraphCut('close', gch);  % do not forget to de-allocate

To see the output labels, you need to reshape

output = reshape(L, size(slices));

PS,
If your spatial distance between slices is larger than the gap between neighboring voxels in the same slice, you might need to use different sig2 for ii and jj that are in the same slice and for ii and jj that are on different slices. This requires a bit of an effort.

Shai
  • 111,146
  • 38
  • 238
  • 371
  • Thanks, I loaded 4 maps into `prob` by for loop , `prob(:,:,i,:)= cat(4, squeeze(pMap.a(1,:,:)), squeeze(pMap.a(2,:,:)), squeeze(pMap.a(3,:,:)), squeeze(pMap.a(4,:,:)))`, however, I am getting error: `Error using GraphCut>OpenGraph (line 383) Wrong size of Dc or Sc`. The size of Dc is 3932160-by-4 after reshaping. I do not know why?! – S.EB Oct 25 '17 at 09:59
  • it's been a while since I last used this package. it might be the case that you need to transpose Dc – Shai Oct 25 '17 at 10:51
  • Thanks a lot, it works, Dc was changed to Dc' however, does not improve the segmentation results that much. I set `sig=0.01`; – S.EB Oct 25 '17 at 12:22
  • you might want to increase the relative weight of the smoothness term. You can do that by multiplying `W` by `\lambda > 1` factor. – Shai Oct 25 '17 at 12:32
  • Thank you very much for your help, I did that one too, however, not that much improvement. Edges became a little bit more smoother. CNN was under-segmenting and I was thinking to improve it by increasing the posterior probability of probability maps by using MRF. – S.EB Oct 25 '17 at 13:03
  • you might need to increase `\lambda` significantly. Try very large value until you see only large smooth blobs in the output, and then decrease to get the proper amount of smoothing you wish to get. – Shai Oct 25 '17 at 13:25
  • Dear Shai, In the demo [http://vision.ucla.edu/~brian/gcmex.html] that Brian Fulkerson provides, the function is taking `CLASS` labels as well `[LABELS ENERGY ENERGYAFTER] = GCMEX(CLASS, UNARY, PAIRWISE, LABELCOST,EXPANSION)` . How it works with your version? Thanks – S.EB Nov 10 '17 at 03:34
  • Could I ask if it is `[gch] = GraphCut('set', gch, labels)`, which `labels` is vectorized ground truth matrix? What is the advantages of giving the labels to GC? Thanks – S.EB Nov 10 '17 at 04:41
  • @S.EB a good initialization can sometimes help in getting a better solution – Shai Nov 10 '17 at 04:55
  • Thanks for your help. May I ask another question, is there any formula to find out `lambda` value? or we should continuously try different values? what do you suggest? Thanks – S.EB Nov 10 '17 at 07:13
  • @S.EB it's an open question. You may find the paper on parametric maxflow by C. Rother et al. relevant – Shai Nov 10 '17 at 07:17