笔记（二）maxflow

1、maxflow.cpp中maxflow3d的python封装

#include <Python.h>
#include "numpy/arrayobject.h"
#include "maxflow-v3.0/graph.h"
#include "util.h"
#include <iostream>
using namespace std;

// example to use numpy object: http://blog.debao.me/2013/04/my-first-c-extension-to-numpy/
// write a c extension ot Numpy: http://folk.uio.no/hpl/scripting/doc/python/NumPy/Numeric/numpy-13.html
static PyObject *
maxflow3d_wrapper(PyObject *self, PyObject *args)
{
    PyObject *I=NULL, *P=NULL, *param=NULL;
    PyArrayObject *arr_I=NULL, *arr_P=NULL;
    if (!PyArg_ParseTuple(args, "OOO", &I, &P, &param)) return NULL;
    
    arr_I = (PyArrayObject*)PyArray_FROM_OTF(I, NPY_FLOAT32, NPY_IN_ARRAY);
    arr_P = (PyArrayObject*)PyArray_FROM_OTF(P, NPY_FLOAT32, NPY_IN_ARRAY);
    if (arr_I == NULL || arr_P == NULL) return NULL;
    float lambda = PyFloat_AsDouble(PyTuple_GET_ITEM(param, 0));
    float sigma  = PyFloat_AsDouble(PyTuple_GET_ITEM(param, 1));
    
    /*vv* code that makes use of arguments *vv*/
    int dimI = PyArray_NDIM(arr_I);             // number of dimensions
    int dimP = PyArray_NDIM(arr_P); 
    npy_intp * shapeI = PyArray_DIMS(arr_I);    // npy_intp array of length nd showing length in each dim
    npy_intp * shapeP = PyArray_DIMS(arr_P);

    if(dimI !=3 && dimI != 4){
        cout << "the input dimension can only be 3 or 4"<<endl;
        return NULL;
    }
    if(dimP != 4){
        cout << "dimension of probabilily map should be 4"<<endl;
        return NULL;
    }
    if(shapeI[0] != shapeP[0] || shapeI[1] != shapeP[1] || shapeI[2] != shapeP[2]){
        cout << "image and probability map have different sizes"<<endl;
        return NULL;
    }
    if(shapeP[3] != 2){
        cout << "probabilily map should have two channels"<<endl;
        return NULL;
    }
    
    int chns = 1;
    if(dimI == 4) chns = shapeI[3];
    npy_intp outshape[3] = {shapeI[0], shapeI[1], shapeI[2]};
    PyArrayObject * arr_L = (PyArrayObject*)  PyArray_SimpleNew(3, outshape, NPY_INT8);
    maxflow3d_inference((unsigned char *) arr_L->data, (const float *) arr_I->data, 
        (const float *) arr_P->data, NULL,
        shapeI[0], shapeI[1], shapeI[2], chns, 2, lambda, sigma);

    Py_DECREF(arr_I);
    Py_DECREF(arr_P);
    Py_INCREF(arr_L);
    return PyArray_Return(arr_L);
}

Matlab and Python wrap of Conditional Random Field (CRF) and fully connected (dense) CRF for 2D and 3D image segmentation, according to the following papers:

[1] Yuri Boykov and Vladimir Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision”, IEEE TPAMI, 2004.

[2] Philipp Krähenbühl and Vladlen Koltun, “Efficient inference in fully connected crfs with gaussian edge potentials”, in NIPS, 2011.

[3] Kamnitsas et al in “Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation”, Medical Image Analysis, 2017.

Dependency
This repository depends on the following packages: Maxflow , DenceCRF and 3D Dense CRF

参考：HiLab-git/SimpleCRF

2、文件引用关系

SimpleCRF模块中封装了maxflow，maxflow.maxflow3d在maxflow.cpp中定义。

maxflow.cpp中

maxflow3d_wrapper(PyObject *self, PyObject *args)

调用

maxflow3d_inference((unsigned char *) arr_L->data, (const float *) arr_I->data, 
        (const float *) arr_P->data, NULL,
        shapeI[0], shapeI[1], shapeI[2], chns, 2, lambda, sigma); 

来自于util.h中定义的：

#include <vector>
#include "maxflow-v3.0/graph.h"
using namespace std;

void maxflow3d_inference(unsigned char * label, const float* img, const float * prob, const unsigned char * seed,
    int D, int H, int W, int chns, int cls, float lambda, float sigma);

util.cpp中实现：

void maxflow3d_inference(unsigned char * label, const float* img, const float * prob, const unsigned char * seed,
    int D, int H, int W, int chns, int cls, float lambda, float sigma)
{
    // currently, only cls == 2 is supported
    typedef Graph<float, float, float> GraphType;
    /*estimated # of nodes*/ /*estimated # of edges*/
    GraphType * g = new GraphType(D*H*W, 2*D*H*W);
    g->add_node(D*H*W);
    float max_weight = -100000;
    for(int x=0; x<D; x++)
    {
        for(int y=0; y<H; y++)
        {
            for(int z=0; z<W; z++)
            {
                std::vector<float> pValue = get_pixel_vector(img, D, H, W, chns, x, y, z);
                std::vector<float> qValue;
                float l2dis, n_weight;
                int pIndex  = x*H*W + y*W + z;
                int qIndex;
                int xn, yn, zn;
                int Xoff[3] = {-1, 0, 0};
                int Yoff[3] = {0, -1, 0};
                int Zoff[3] = {0, 0, -1};
                for(int i=0; i<3; i++){
                    xn = x + Xoff[i];
                    yn = y + Yoff[i];
                    zn = z + Zoff[i];
                    if(xn < 0 || yn < 0 || zn < 0) continue;
                    qValue   = get_pixel_vector(img, D, H, W, chns, xn, yn, zn);
                    l2dis    = get_l2_distance(pValue, qValue);
                    n_weight = lambda*exp(-l2dis * l2dis/(2*sigma*sigma));
                    qIndex   = xn*H*W + yn*W + zn;
                    g->add_edge(qIndex,pIndex,n_weight,n_weight);
                    if(n_weight > max_weight) max_weight = n_weight;
                }
            }
        }
    }

    max_weight =  1000 * max_weight;
    for(int x=0; x<D; x++)
    {
        for(int y=0; y<H; y++)
        {
            for(int z=0; z<W; z++)
            {
                bool seed_exist = false;
                float s_weight  = 1e-3;
                float t_weight  = 1e-3;
                if(seed != NULL)
                {
                    std::vector<unsigned char> label = get_pixel_vector(seed, D, H, W, 2, x, y, z); 
                    if(label[0] > 0){
                        t_weight   = max_weight;
                        seed_exist = true; 
                    }
                    else if(label[1] > 0){
                        s_weight   = max_weight;
                        seed_exist = true;
                    }
                }
                if(!seed_exist){
                    std::vector<float> probs = get_pixel_vector(prob, D, H, W, 2, x, y, z);
                    s_weight = -log(probs[0]);
                    t_weight = -log(probs[1]);
                }
                int  pIndex  = x*H*W + y*W + z;
                g->add_tweights(pIndex,s_weight,t_weight);
            }
        }
    }
    double flow = g->maxflow();
    // cout<<"max flow: "<<flow<<endl;
    for(int i=0; i<D*H*W; i++) label[i] = 1 - g->what_segment(i);
    delete g;
}

maxflow是图割算法，它可以被GraphCut调用，也可以被GrabCut调用。

参考：
参考：Source code for monailabel.scribbles.utils