2014年3月29日 星期六

I've been awarded research grant: 748,000SGD, for research topic: Quantitative Tumour Angiogenesis Assessments using 3D Radiographic/Histopathologic Imaging and Xenotransplantation for Cancer Therapy...

張貼者： 於 沒有留言:

2012年2月18日 星期六

TV-L1 based Optical Flow using GPU

Mr. Chris Mcclanahan has implemented a TV-L1 based Optical Flow based on LibJacket (see here). However, is is not compatible the newer version of LibJacket...well...there is no more LibJacket. Now AccelerEyes announced ArrayFire, which has a huge gap to LibJacket. Thus, I managed modified his work as the following:

//
// Chris McClanahan - 2011
// Modified by Chao-Hui Huang 2012
//
// Adapted from: http://gpu4vision.icg.tugraz.at/index.php?content=downloads.php
//   "An Improved Algorithm for TV-L1 Optical Flow"
//
// More info: http://mcclanahoochie.com/blog/portfolio/gpu-tv-l1-optical-flow-with-libjacket/
//

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <string.h>

#include <opencv/cv.h>
#include <opencv/cxcore.h>
#include <opencv/highgui.h>

#include <arrayfire.h>

using namespace std;
using namespace cv;
using namespace af;

// control
const float pfactor = 0.7;    // scale each pyr level by this amount
const int max_plevels = 9;    // number of pyramid levels
const int max_iters = 6;      // u v w update loop
const float lambda = 40;      // smoothness constraint
const int max_warps = 3;      // warping u v warping
const int min_img_sz = 20;    // min mxn img in pyramid
#define TIMING 0              // warmup, then average multiple runs

// functions
int  grab_frame(Mat& img, char* filename);
void create_pyramids(array& im1, array& im2, array& pyr1, array& pyr2);
void process_pyramids(array& pyr1, array& pyr2, array& u, array& v);
void tv_l1_dual(array& u, array& v, array& p, array& w, array& I1, array& I2, int level);
void optical_flow_tvl1(Mat& img1, Mat& img2, Mat& u, Mat& v);
void display_flow(array& I2, array& u, array& v);
void MatToFloat(const Mat& thing, float* thing2);
void FloatToMat(float const* thing, Mat& thing2);

// misc
int plevels = max_plevels;
const int n_dual_vars = 6;
static int cam_init = 0;
static int pyr_init = 0;
VideoCapture  capture;
int pyr_M[max_plevels + 1];
int pyr_N[max_plevels + 1];
array pyr1, pyr2;

// macros
#define MSG(msg,...) do {                                   \
                fprintf(stdout,__FILE__":%d(%s) " msg "\n",     \
                        __LINE__, __FUNCTION__, ##__VA_ARGS__); \
                fflush(stdout);                                 \
            } while (0)

#define M_PI 3.14159265358979323846

// ===== main =====
void optical_flow_tvl1(Mat& img1, Mat& img2, Mat& mu, Mat& mv) {

    // extract cv image 1
    Mat mi1(img1.rows, img1.cols, CV_8UC1);
    cvtColor(img1.t(), mi1, CV_BGR2GRAY);
    mi1.convertTo(mi1, CV_32FC1);
    float* fi1 = (float*)mi1.data;
    array I1 = array(img1.rows, img1.cols, fi1) / 255.0f;

    // extract cv image 2
    Mat mi2(img2.rows, img2.cols, CV_8UC1);
    cvtColor(img2.t(), mi2, CV_BGR2GRAY);
    mi2.convertTo(mi2, CV_32FC1);
    float* fi2 = (float*)mi2.data;
    array I2 = array(img2.rows, img2.cols, fi2) / 255.0f;

#if TIMING
    // runs
    int nruns = 4;
    // warmup
    create_pyramids(I1, I2, pyr1, pyr2);
    f32 ou, ov;
    process_pyramids(pyr1, pyr2, ou, ov);
    // timing
    timer::tic();
    for (int i = 0; i < nruns; ++i) {
        create_pyramids(I1, I2, pyr1, pyr2);
        process_pyramids(pyr1, pyr2, ou, ov);
    }
    MSG("fps: %f", 1.0f / (timer::toc() / (float)nruns));
#else
    // timing
    timer::tic();
    // pyramids
    create_pyramids(I1, I2, pyr1, pyr2);
    // flow
    array ou, ov;
    process_pyramids(pyr1, pyr2, ou, ov);
    // timing
    MSG("fps: %f", 1.0f / (timer::toc()));
#endif

    // output
#if 1
    // to opencv
    FloatToMat((float*)ou.T().host(), mu);
    FloatToMat((float*)ov.T().host(), mv);
#else
    // to libjacket
    display_flow(I2, ou, ov);
#endif
}


void MatToFloat(const Mat& thing, float* thing2) {
    int tmp = 0;
    for (int i = 0; i < thing.rows; i++) {
        const float* fptr = thing.ptr<float>(i);
        for (int j = 0; j < thing.cols; j++)
        { thing2[tmp++] = fptr[j]; }
    }
}


void FloatToMat(float const* thing, Mat& thing2) {
    int tmp = 0;
    for (int i = 0; i < thing2.rows; ++i) {
        float* fptr = thing2.ptr<float>(i);
        for (int j = 0; j < thing2.cols; ++j)
        { fptr[j] = thing[tmp++]; }
    }
}


void display_flow(array& I2, array& u, array& v) {
#if 1
    // show in libjacket
    palette("bone");
    subfigure(2, 2, 1); imgplot(I2);                  title("input");
    subfigure(2, 2, 2); imgplot(u);                   title("u");
    subfigure(2, 2, 3); imgplot(v);                   title("v");
    subfigure(2, 2, 4); imgplot((abs(v) + abs(u)));   title("u+v");
    // int M = I2.dims()[0];
    // int N = I2.dims()[1];
    // f32 idx, idy; meshgrid(idx, idy, f32(seq(0,N-1,3)), f32(seq(0,M-1,3)));
    // quiver(idx,idy,u,v);
    draw();
#else
    // show in opencv
    int M = I2.dims()[0];
    int N = I2.dims()[1];
    Mat mu(M, N, CV_32FC1);
    Mat mv(M, N, CV_32FC1);
    FloatToMat(u.T().host(), mu);
    FloatToMat(v.T().host(), mv);
    imshow("u", mu);
    imshow("v", mv);
#endif
}


void display_flow(const Mat& u, const Mat& v) {
#if 0
    cv::Mat magnitude, angle, bgr;
    cv::cartToPolar(u, v, magnitude, angle, true);
    double mag_max, mag_min;
    cv::minMaxLoc(magnitude, &mag_min, &mag_max);
    magnitude.convertTo(magnitude, -1, 1.0 / mag_max);
    cv::Mat _hsv[3], hsv_image;
    _hsv[0] = angle;
    _hsv[1] = Mat::ones(angle.size(), CV_32F);
    _hsv[2] = magnitude;
    cv::merge(_hsv, 3, hsv_image);
#else
    cv::Mat magnitude, angle, bgr;
    Mat hsv_image(u.rows, u.cols, CV_8UC3);
    for (int i = 0; i < u.rows; ++i) {
        const float* x_ptr = u.ptr<float>(i);
        const float* y_ptr = v.ptr<float>(i);
        uchar* hsv_ptr = hsv_image.ptr<uchar>(i);
        for (int j = 0; j < u.cols; ++j, hsv_ptr += 3, ++x_ptr, ++y_ptr) {
            hsv_ptr[0] = (uchar)((atan2f(*y_ptr, *x_ptr) / M_PI + 1) * 90);
            hsv_ptr[1] = hsv_ptr[2] = (uchar) std::min<float>(
                                          sqrtf(*y_ptr * *y_ptr + *x_ptr * *x_ptr) * 20, 255.0);
        }
    }
#endif
    cv::cvtColor(hsv_image, bgr, CV_HSV2BGR);
    cv::imshow("optical flow", bgr);
}


int grab_frame(Mat& img, char* filename) {

    // camera/image setup
    if (!cam_init) {
        if (filename != NULL) {
            capture.open(filename);
        } else {
            float rescale = 0.615;
            int w = 640 * rescale;
            int h = 480 * rescale;
            capture.open(0); //try to open
            capture.set(CV_CAP_PROP_FRAME_WIDTH, w);  capture.set(CV_CAP_PROP_FRAME_HEIGHT, h);
        }
        if (!capture.isOpened()) { cerr << "open video device fail\n" << endl; return 0; }
        capture >> img; capture >> img;
        if (img.empty()) { cout << "load image fail " << endl; return 0; }
        namedWindow("cam", CV_WINDOW_KEEPRATIO);
        printf(" img = %d x %d \n", img.cols, img.rows);
        cam_init = 1;
    }

    // get frames
    capture.grab();
    capture.retrieve(img);
    imshow("cam", img);

    if (waitKey(10) >= 0) { return 0; }
    else { return 1; }
}


void gen_pyramid_sizes(array& im1) {
    dim4 mnk = im1.dims();
    float sM = mnk[0];
    float sN = mnk[1];
    // store resizing
    for (int level = 0; level <= plevels; ++level) {
        if (level == 0) {
        } else {
            sM *= pfactor;
            sN *= pfactor;
        }
        pyr_M[level] = (int)(sM + 0.5f);
        pyr_N[level] = (int)(sN + 0.5f);
        MSG(" pyr %d: %d x %d ", level, (int)sM, (int)sN);
        if (sM < min_img_sz || sN < min_img_sz) { plevels = level; break; }
    }
}

void create_pyramids(array& im1, array& im2, array& pyr1, array& pyr2) {

    if (!pyr_init) {
        // list of h,w
        gen_pyramid_sizes(im1);

        // init
        pyr1 = zeros(pyr_M[0], pyr_N[0], plevels);
        pyr2 = zeros(pyr_M[0], pyr_N[0], plevels);
        pyr_init = 1;
    }

    // create
    for (int level = 0; level < plevels; level++) {
        if (level == 0) {
            pyr1(span, span, level) = im1;
            pyr2(span, span, level) = im2;
        } else {
            seq spyi = seq(pyr_M[level - 1]);
            seq spxi = seq(pyr_N[level - 1]);
            array small1 = resize(pyr1(spyi, spxi, level - 1), pyr_M[level], pyr_N[level], AF_RSZ_Bilinear);
            array small2 = resize(pyr2(spyi, spxi, level - 1), pyr_M[level], pyr_N[level], AF_RSZ_Bilinear);
            seq spyo = seq(pyr_M[level]);
            seq spxo = seq(pyr_N[level]);
            pyr1(spyo, spxo, level) = small1;
            pyr2(spyo, spxo, level) = small2;
        }
    }
}


void process_pyramids(array& pyr1, array& pyr2, array& ou, array& ov) {
    array p, u, v, w;

    // pyramid loop
    for (int level = plevels - 1; level >= 0; level--) {
        if (level == plevels - 1) {
            u  = zeros(pyr_M[level], pyr_N[level]);
            v  = zeros(pyr_M[level], pyr_N[level]);
            w  = zeros(pyr_M[level], pyr_N[level]);
            p  = zeros(pyr_M[level], pyr_N[level], n_dual_vars);
        } else {
            float rescale_u =  pyr_N[level + 1] / (float)pyr_N[level];
            float rescale_v =  pyr_M[level + 1] / (float)pyr_M[level];
            // propagate
            array u_ =  resize(u, pyr_M[level], pyr_N[level], AF_RSZ_Bilinear) * rescale_u;
            array v_ =  resize(v, pyr_M[level], pyr_N[level], AF_RSZ_Bilinear) * rescale_v;
            array w_ =  resize(w, pyr_M[level], pyr_N[level], AF_RSZ_Bilinear);
            array p_ = zeros(pyr_M[level], pyr_N[level], n_dual_vars);
            gfor(array ndv, n_dual_vars) {
                p_(span, span, ndv) = resize(p(span, span, ndv), pyr_M[level], pyr_N[level], AF_RSZ_Bilinear);
            }
            u = u_;  v = v_;  p = p_;  w = w_;
        }

        // extract
        seq spy = seq(pyr_M[level]);
        seq spx = seq(pyr_N[level]);
        array I1 = pyr1(spy, spx, level);
        array I2 = pyr2(spy, spx, level);

        // ===== core ====== //
        tv_l1_dual(u, v, p, w, I1, I2, level);
        // ===== ==== ====== //
    }

    // output
    ou = u;
    ov = v;
}


void warping(array& Ix, array& Iy, array& It, array& I1, array& I2, array& u, array& v) {

    dim4 mnk = I2.dims();
    int M = mnk[0];
    int N = mnk[1];
    array idx = tile(array(seq(N)).T(), M, 1) + 1;
    array idy = tile(array(seq(M)), 1, N) + 1;
    /* ^ BUG: idx idy should ideally be [0-N); ^ */

    array idxx0 = idx + u;
    array idyy0 = idy + v;
    array idxx = max(1, min(N - 1, idxx0));
    array idyy = max(1, min(M - 1, idyy0));

    // interp2 based warp ()
    It = interp(idy, idx, I2, idyy, idxx) - I1;

    // interp2 based warp ()
    array idxm = max(1, min(N - 1, idxx - 1.f));
    array idxp = max(1, min(N - 1, idxx + 1.f));
    array idym = max(1, min(M - 1, idyy - 1.f));
    array idyp = max(1, min(M - 1, idyy + 1.f));
    Ix = interp(idy, idx, I2, idy, idxp) - interp(idy, idx, I2, idy, idxm);
    Iy = interp(idy, idx, I2, idyp, idx) - interp(idy, idx, I2, idym, idx);
    /* ^ BUG: interp2 should be cubic; that may fix things; ^ */
}


void dxym(array& Id, array I0x, array I0y) {
    // divergence
    dim4 mnk = I0x.dims();
    int M = mnk[0];
    int N = mnk[1];

    array x0 = zeros(M, N);
    array x1 = zeros(M, N);
    x0(seq(M - 1), seq(N)) = I0x(seq(M - 1), seq(N));
    x1(seq(1,M-1), seq(N)) = I0x(seq(1,M-1), seq(N));

    array y0 = zeros(M, N);
    array y1 = zeros(M, N);
    y0(seq(M), seq(N - 1)) = I0y(seq(M), seq(N - 1));
    y1(seq(M), seq(1,N-1)) = I0y(seq(M), seq(1,N-1));

    Id = (x0 - x1) + (y0 - y1);
}


void dxyp(array& Ix, array& Iy, array& I0) {
    // shifts
    dim4 mnk = I0.dims();
    int M = mnk[0];
    int N = mnk[1];

    array y0 = I0;
    array y1 = I0;
    y0(seq(0, M - 2), span) = I0(seq(1, M - 1), span);

    array x0 = I0;
    array x1 = I0;
    x0(span, seq(0, N - 2)) = I0(span, seq(1, N - 1));

    Ix = (x0 - x1);  Iy = (y0 - y1);
}


void tv_l1_dual(array& u, array& v, array& p, array& w, array& I1, array& I2, int level) {
    try {
    float L = sqrtf(8.0f);
    float tau   = 1 / L;
    float sigma = 1 / L;

    float eps_u = 0.01f;
    float eps_w = 0.01f;
    float gamma = 0.02f;

    array u_ = u;
    array v_ = v;
    array w_ = w;

    for (int j = 0; j < max_warps; j++) {

        array u0 = u;
        array v0 = v;

        // warping
        array Ix, Iy, It;   warping(Ix, Iy, It, I1, I2, u0, v0);

        // gradients
        array I_grad_sqr = max(float(1e-6), array(pow(Ix, 2) + pow(Iy, 2) + gamma * gamma));

        // inner loop
        for (int k = 0; k < max_iters; ++k) {

            // dual =====

            // shifts
            array u_x, u_y;    dxyp(u_x, u_y, u_);
            array v_x, v_y;    dxyp(v_x, v_y, v_);
            array w_x, w_y;    dxyp(w_x, w_y, w_);

            // update dual
            p(span, span, 0) = (p(span, span, 0) + sigma * u_x) / (1 + sigma * eps_u);
            p(span, span, 1) = (p(span, span, 1) + sigma * u_y) / (1 + sigma * eps_u);
            p(span, span, 2) = (p(span, span, 2) + sigma * v_x) / (1 + sigma * eps_u);
            p(span, span, 3) = (p(span, span, 3) + sigma * v_y) / (1 + sigma * eps_u);

            p(span, span, 4) = (p(span, span, 4) + sigma * w_x) / (1 + sigma * eps_w);
            p(span, span, 5) = (p(span, span, 5) + sigma * w_y) / (1 + sigma * eps_w);

            // normalize
            array reprojection = max(1, sqrt(pow(p(span, span, 0), 2) + pow(p(span, span, 1), 2) +
                                           pow(p(span, span, 2), 2) + pow(p(span, span, 3), 2)));

            p(span, span, 0) = p(span, span, 0) / reprojection;
            p(span, span, 1) = p(span, span, 1) / reprojection;
            p(span, span, 2) = p(span, span, 2) / reprojection;
            p(span, span, 3) = p(span, span, 3) / reprojection;

            reprojection = max(1, sqrt(pow(p(span, span, 4), 2) + pow(p(span, span, 5), 2)));

            p(span, span, 4) = p(span, span, 4) / reprojection;
            p(span, span, 5) = p(span, span, 5) / reprojection;

            // primal =====

            // divergence
            array div_u;   dxym(div_u, p(span, span, 0), p(span, span, 1));
            array div_v;   dxym(div_v, p(span, span, 2), p(span, span, 3));
            array div_w;   dxym(div_w, p(span, span, 4), p(span, span, 5));

            // old
            u_ = u;
            v_ = v;
            w_ = w;

            // update
            u = u + tau * div_u;
            v = v + tau * div_v;
            w = w + tau * div_w;

            // indexing
            array rho  = It + mul((u - u0), Ix) + mul((v - v0), Iy) + gamma * w;
            array idx1 = rho      <  -tau * lambda * I_grad_sqr;
            array idx2 = rho      >   tau * lambda * I_grad_sqr;
            array idx3 = abs(rho) <=  tau * lambda * I_grad_sqr;

            u = u + tau * lambda * (mul(Ix, idx1)) ;
            v = v + tau * lambda * (mul(Iy, idx1)) ;
            w = w + tau * lambda * gamma * idx1;

            u = u - tau * lambda * (mul(Ix, idx2)) ;
            v = v - tau * lambda * (mul(Iy, idx2)) ;
            w = w - tau * lambda * gamma * idx2;

            u = u - mul(rho, mul(idx3, Ix / I_grad_sqr));
            v = v - mul(rho, mul(idx3, Iy / I_grad_sqr));
            w = w - mul(rho, mul(idx3, gamma / I_grad_sqr));

            // propagate
            u_ = 2 * u - u_;
            v_ = 2 * v - v_;
            w_ = 2 * w - w_;

        }

        // output
        // const unsigned hw[] = {3, 3};
        u = medfilt(u, 3, 3);
        v = medfilt(v, 3, 3);

    } /* j < warps */


    } catch (af::exception& e) {
        cout << e.what() << endl;
        throw;
    }


}


// =======================================

int main(int argc, char* argv[]) {

    // video file or usb camera
    Mat cam_img, prev_img, disp_u, disp_v;
    int is_images = 0;
    if (argc == 2) { grab_frame(prev_img, argv[1]); } // video
    else if (argc == 3) {
        prev_img = imread(argv[1]); cam_img = imread(argv[2]);
        is_images = 1;
    } else { grab_frame(prev_img, NULL); } // usb camera

    // results
    int mm = prev_img.rows;  int nn = prev_img.cols;
    disp_u = Mat::zeros(mm, nn, CV_32FC1);
    disp_v = Mat::zeros(mm, nn, CV_32FC1);
    printf("img %d x %d \n", mm, nn);

    // process main
    if (is_images) {
        // show
        imshow("i", cam_img);
        // process files
        optical_flow_tvl1(prev_img, cam_img, disp_u, disp_v);
        // show
        // imshow("u", disp_u);
        // imshow("v", disp_v);
        display_flow(disp_u, disp_v);
        waitKey(0);
        // // write
        // writeFlo(disp_u, disp_v);
    } else {
        // process loop
        while (grab_frame(cam_img, NULL)) {
            try {
                // process
                optical_flow_tvl1(prev_img, cam_img, disp_u, disp_v);
                // frames
                prev_img = cam_img.clone();
                // show
                // imshow("u", disp_u);
                // imshow("v", disp_v);
                display_flow(disp_u, disp_v);
            } catch (af::exception& e) {
                cout << e.what() << endl;
                throw;
            }
        }
    }

    return 0;
}
張貼者： 於 沒有留言:

2011年3月17日 星期四

The page of research results is updated...

Check it out...

http://web.bii-sg.org/~huangch/significant_results/index.html
張貼者： 於 沒有留言:

2011年3月10日 星期四

Some Results of In-vitro Neural Stem Cell Tracking and Acquisition

I have posted the results of In-vitro Neural Stem Cell Tracking and Acquisition...

http://www.youtube.com/watch?v=Dzj7wBVOPzs

http://www.youtube.com/watch?v=jDDwpp9FYMI

Check it out...
張貼者： 於 沒有留言:

2011年3月1日 星期二

Some good results in our project: Consciousness4UAS

The aim of this project is to discover the visual consciousness model and its applications in the unmanned aircraft system.

We introduced an model for consciousness based visual attention model, and applied on the aerial video (simulating the real time images obtained from an UAV)

Results: http://www.youtube.com/user/consciousness4uas
張貼者： 於 沒有留言:
訂閱： 文章 (Atom)