cd: just add current code, regardless of hygiene.

cd.c

@@ -1,21 +1,35 @@
 // processing changedetection.net dataset
 #include <stdio.h>
 #include <stdlib.h>
+#include <pthread.h>
 
 #include <opencv2/imgproc/imgproc_c.h>
 #include <opencv2/highgui/highgui_c.h>
 
 #include "kdtree.h"
 #include "prop.h"
 #include "wht.h"
+#include "sal.h"
+#include "gck.h"
 
 #define XY_TO_INT(x, y) (((y) << 16) | (x))
 #define XY_TO_X(x) ((x)&((1<<16)-1))
 #define XY_TO_Y(y) ((y)>>16)
 
+#define KERNS 8
+
+#include <sys/time.h>
+static inline double get_time()
+{
+    struct timeval t;
+    gettimeofday(&t, NULL);
+    return t.tv_sec + t.tv_usec * 1e-6;
+}
+
 static void print_usage(char **argv)
 {
-    printf("Usage: %s <path> <start> <end>\n", argv[0]);
+    printf("Usage: %s <path> <start> <end> <bkg_name> <diff_name>\n",
+        argv[0]);
     exit(1);
 }
 
@@ -139,21 +153,21 @@ static int* block_coeffs(IplImage *img, int* plane_coeffs) {
     int dim = plane_coeffs[0] + plane_coeffs[1] + plane_coeffs[2];
     int sz = size.width*size.height/64*dim;
     int *buf = malloc(sizeof(int)*sz);
-    unsigned *order_p0 = build_path(plane_coeffs[0], 8);
-    unsigned *order_p1 = build_path(plane_coeffs[1], 8);
-    unsigned *order_p2 = build_path(plane_coeffs[2], 8);
+    unsigned *order_p0 = build_path(plane_coeffs[0], KERNS);
+    unsigned *order_p1 = build_path(plane_coeffs[1], KERNS);
+    unsigned *order_p2 = build_path(plane_coeffs[2], KERNS);
 
     cvSplit(img, b, g, r, NULL);
 
     wht2d(b, trans);
-    quantize(trans, plane_coeffs[0], 8, order_p0, buf, dim);
+    quantize(trans, plane_coeffs[0], KERNS, order_p0, buf, dim);
 
     wht2d(g, trans);
-    quantize(trans, plane_coeffs[1], 8, order_p1,
+    quantize(trans, plane_coeffs[1], KERNS, order_p1,
         buf+plane_coeffs[0], dim);
 
     wht2d(r, trans);
-    quantize(trans, plane_coeffs[2], 8, order_p2,
+    quantize(trans, plane_coeffs[2], KERNS, order_p2,
         buf+plane_coeffs[0]+plane_coeffs[1], dim);
 
     cvReleaseImage(&trans);
@@ -176,19 +190,19 @@ static IplImage* splat(int *coeffs, CvSize size, int *plane_coeffs)
     IplImage *img = cvCreateImage(size, IPL_DEPTH_8U, 3);
     IplImage *trans = cvCreateImage(size, IPL_DEPTH_16S, 1);
     int dim = plane_coeffs[0] + plane_coeffs[1] + plane_coeffs[2];
-    unsigned *order_p0 = build_path(plane_coeffs[0], 8);
-    unsigned *order_p1 = build_path(plane_coeffs[1], 8);
-    unsigned *order_p2 = build_path(plane_coeffs[2], 8);
+    unsigned *order_p0 = build_path(plane_coeffs[0], KERNS);
+    unsigned *order_p1 = build_path(plane_coeffs[1], KERNS);
+    unsigned *order_p2 = build_path(plane_coeffs[2], KERNS);
 
     memset(trans->imageData, 0, trans->imageSize);
-    dequantize(trans, plane_coeffs[0], order_p0, 8, coeffs, dim);
+    dequantize(trans, plane_coeffs[0], order_p0, KERNS, coeffs, dim);
     iwht2d(trans, g);
     memset(trans->imageData, 0, trans->imageSize);
-    dequantize(trans, plane_coeffs[1], order_p1, 8,
+    dequantize(trans, plane_coeffs[1], order_p1, KERNS,
         coeffs+plane_coeffs[0], dim);
     iwht2d(trans, b);
     memset(trans->imageData, 0, trans->imageSize);
-    dequantize(trans, plane_coeffs[2], order_p2, 8,
+    dequantize(trans, plane_coeffs[2], order_p2, KERNS,
         coeffs+plane_coeffs[0]+plane_coeffs[1], dim);
     iwht2d(trans, r);
 
@@ -263,13 +277,16 @@ static IplImage *alignedImage(CvSize dim, int depth, int chan, int align)
     return cvCreateImage(s, depth, chan);
 }
 
+int orig_w = -1, orig_h = -1;
 static IplImage *alignedImageFrom(char *file, int align)
 {
     IplImage *pre = cvLoadImage(file, CV_LOAD_IMAGE_COLOR);
     IplImage *img = alignedImage(cvGetSize(pre), pre->depth, pre->nChannels, align);
     char *pre_data = pre->imageData;
     char *img_data = img->imageData;
     int i;
+    orig_w = pre->width;
+    orig_h = pre->height;
     for (i = 0; i < pre->height; i++) {
         memcpy(img_data, pre_data, pre->widthStep);
         img_data += img->widthStep;
@@ -280,47 +297,246 @@ static IplImage *alignedImageFrom(char *file, int align)
 }
 
 static int plane_coeffs[] = {2, 9, 5};
-IplImage *recon_g, *diff_g;
-static void process(kd_tree *kdt, IplImage *bkg, IplImage *img)
+IplImage *recon_g, *diff_g, *lap_g, *gray_g, *mask_g, *bkg_g;
+static void init_g(IplImage *bkg)
+{
+    CvSize bsz = cvGetSize(bkg);
+    recon_g = cvCreateImage(bsz, bkg->depth, bkg->nChannels);
+    diff_g  = cvCreateImage(bsz, bkg->depth, bkg->nChannels);
+    bkg_g   = cvCreateImage(bsz, bkg->depth, bkg->nChannels);
+    gray_g  = cvCreateImage(bsz, bkg->depth, 1);
+    mask_g  = cvCreateImage(bsz, bkg->depth, 1);
+    lap_g   = cvCreateImage(bsz, IPL_DEPTH_32F, bkg->nChannels);
+}
+
+static void free_g()
 {
-    int *imgc = block_coeffs(img, plane_coeffs);
-    CvSize img_blks = {(img->width/8)+7, (img->height/8)+7};
-    IplImage *xy = prop_match_complete(kdt, imgc, bkg, img_blks);
+    cvReleaseImage(&recon_g);
+    cvReleaseImage(&diff_g);
+    cvReleaseImage(&gray_g);
+    cvReleaseImage(&mask_g);
+    cvReleaseImage(&lap_g);
+    cvReleaseImage(&bkg_g);
+}
+
+static void process(kd_tree *kdt, IplImage *bkg, IplImage *diff,
+    IplImage *img, char *outname)
+{
+    //int *imgc = block_coeffs(img, plane_coeffs);
+    cvAbsDiff(img, bkg, diff_g);
+    int *imgc = block_coeffs(diff_g, plane_coeffs);
+    CvSize blksz = {(img->width/8)+7, (img->height/8)+7};
+    IplImage *xy = prop_match_complete(kdt, imgc, bkg, blksz);
     IplImage *rev = splat(imgc, cvGetSize(img), plane_coeffs);
-    xy2blks_special(xy, bkg, recon_g, 8);
-    cvAbsDiff(recon_g, rev, diff_g);
+    xy2blks_special(xy, diff, recon_g, 8);
+    cvAbsDiff(rev, recon_g, diff_g);
+    cvReleaseImage(&rev);
+    /*int *imgc;
+    prop_coeffs(diff_g, plane_coeffs, &imgc);
+    CvSize blksz = cvGetSize(bkg);
+    IplImage *xy = prop_match_complete(kdt, imgc, bkg, blksz);
+    xy2img(xy, diff, recon_g);
+    cvAbsDiff(diff_g, recon_g, diff_g);*/
+    //cvShowImage("diff_g before mul", diff_g);
+    //cvAbsDiff(diff_g, diff, diff_g);
+    //cvMul(diff_g, idiff, diff_g, 1);
+    cvCvtColor(diff_g, gray_g, CV_BGR2GRAY);
+
+    // full pixel dc
+    //IplImage *dc = make_dc(gray_g);
+    IplImage *mask = cvCreateImage(cvGetSize(img), IPL_DEPTH_8U, 1);
+    IplImage *dc_f = cvCreateImage(cvGetSize(img), IPL_DEPTH_32F, 1);
+    int *graydc = gck_calc_2d((uint8_t*)gray_g->imageData, img->width, img->height, KERNS, 1);
+    IplImage *dc = cvCreateImageHeader(cvGetSize(img), IPL_DEPTH_32S, 1);
+    int step = img->width + KERNS - 1;
+    cvSetData(dc, graydc, step*sizeof(int));
+    cvConvertScale(dc, dc_f, 1/255.0, 0);
+    cvThreshold(gray_g, mask, 25, 255.0, CV_THRESH_BINARY);
+
+    mask->width = orig_w;
+    mask->height = orig_h;
+    if (*outname) cvSaveImage(outname, mask, 0);
+
+    /*double min = 0, max = 0;
+    cvMinMaxLoc(dc, &min, &max, NULL, NULL, NULL);
+    printf("min: %3f max: %3f\n", min, max);
+    CvScalar scalar = cvRealScalar(-min);
+    cvAddS(dc, scalar, dc_f, NULL);
+    cvConvertScale(dc_f, dc_f, 1.0/(max - min), 0);*/
+
+    // macroblock based counts
+    //cvAdaptiveThreshold(gray_g, mask_g, 255, CV_ADAPTIVE_THRESH_GAUSSIAN_C, CV_THRESH_BINARY, 3, 0);
+
+    //cvSmooth(mask_g, mask_g, CV_GAUSSIAN, 9, 9, 0, 0);
+    //cvSmooth(diff_g, diff_g, CV_GAUSSIAN, 5, 5, 0, 0);
+    //cvLaplace(diff_g, lap_g, 3);
     //IplImage *xy = prop_match(bkg, img);
     //xy2img(xy, bkg, recon_g);
     //cvAbsDiff(recon_g, img, diff_g);
-    cvShowImage("img", recon_g);
-    cvShowImage("diff", diff_g);
+    /*cvShowImage("recon", recon_g);
+    //cvShowImage("diff", rev);
+    cvShowImage("diff_g", diff_g);
+    cvShowImage("img", img);
+    cvShowImage("gray", gray_g);
+    cvShowImage("dc float", dc_f);
+    cvShowImage("mask", mask);
+    cvShowImage("dc", dc);*/
     free(imgc);
     cvReleaseImage(&xy);
+    cvReleaseImageHeader(&dc);
+    free(graydc);
+    cvReleaseImage(&mask);
+    cvReleaseImage(&dc_f);
+}
+
+CvHistogram *make_hist(IplImage *img)
+{
+    int numBins = 256;
+    float range[] = {0, 255};
+    float *ranges[] = { range };
+    CvHistogram *hist = cvCreateHist(1, &numBins, CV_HIST_ARRAY, ranges, 1);
+    cvCalcHist(&img, hist, 0, 0);
+    return hist;
+}
+
+static void test(IplImage *b)
+{
+    CvSize s = cvGetSize(b);
+    //s.width /= 2;
+    //s.height /= 2;
+    IplImage *rsz = alignedImage(s, b->depth, b->nChannels, 8);
+    cvResize(b, rsz, CV_INTER_CUBIC);
+    IplImage *sal = saliency(rsz);
+    cvShowImage("saliency", sal);
+    cvReleaseImage(&rsz);
+    cvReleaseImage(&sal);
+    return;
+}
+
+typedef struct thread_ctx {
+    int start;
+    int end;
+    int nb;
+    char *path;
+    char *outfile;
+    IplImage *bkg;
+    IplImage *diff;
+    kd_tree *kdt;
+} thread_ctx;
+
+static void* run_thr(void *arg)
+{
+    thread_ctx *ctx = (thread_ctx*)arg;
+    double t = 0;
+    int i, start = ctx->start, end = ctx->end, nb = ctx->nb;
+    char outname[1024], *path = ctx->path, *outfile = ctx->outfile;
+    memset(outname, '\0', sizeof(outname));
+    printf("starting thread %d\n", start);
+    for (i = start; i <= end; i+= nb) {
+        IplImage *img = alignedImageFrom(mkname(path, i), 8);
+        double start = get_time();
+        if (outfile) snprintf(outname, sizeof(outname), "%s/bin%06d.png", outfile, i);
+        process(ctx->kdt, ctx->bkg, ctx->diff, img, outname);
+        t += (get_time() - start);
+        if ((cvWaitKey(1)&255)==27)break; // esc
+        cvReleaseImage(&img);
+    }
+    printf("avg: %fms\n", t/(i-start)*1000);
+    return NULL;
 }
 
 int main(int argc, char **argv)
 {
-    if (argc < 4) print_usage(argv);
+    if (argc < 6) print_usage(argv);
     char *path = argv[1];
     int start = atoi(argv[2]), end = atoi(argv[3]), i, *bkgc;
-    IplImage *bkg = alignedImageFrom(mkname(path, start-1), 8);
+    //IplImage *bkg = alignedImageFrom(mkname(path, 1), 8);
+    IplImage *bkg = alignedImageFrom(argv[4], 8);
     int dim = plane_coeffs[0] + plane_coeffs[1] + plane_coeffs[2];
     CvSize bsz = cvGetSize(bkg);
+    IplImage *d8 = alignedImageFrom(argv[5], 8);
+    //IplImage *d8 = alignedImageFrom(mkname(path, 1), 8);
     int w = bsz.width - 8 + 1, h = bsz.height - 8 + 1, sz = w*h;
     kd_tree kdt;
 
-    printf("cd: %s %d %d\n", path, start, end);
+    printf("cd: %s %d %d %s %s\n", path, start, end, argv[4], argv[5]);
+    init_g(bkg);
     memset(&kdt, 0, sizeof(kd_tree));
+
+    /*
+    IplImage *b32 = cvCreateImage(bsz, IPL_DEPTH_32F, bkg->nChannels);
+    IplImage *i32 = cvCreateImage(bsz, IPL_DEPTH_32F, bkg->nChannels);
+    IplImage *d32 = cvCreateImage(bsz, IPL_DEPTH_32F, bkg->nChannels);
+    IplImage *diff= cvCreateImage(bsz, IPL_DEPTH_32F, bkg->nChannels);
+    cvXor(b32, b32, b32, NULL);
+    cvXor(d32, d32, d32, NULL);
+    cvConvertScale(bkg, b32, 1/255.0, 0);
+    for (i = 1; i < start; i++) {
+        IplImage *img = alignedImageFrom(mkname(path, i), 8);
+        cvConvertScale(img, i32, 1/256.0, 0);
+        cvAbsDiff(i32, b32, diff);
+        cvRunningAvg(diff, d32, 1.0/start, NULL);
+        cvRunningAvg(i32, b32, 1.0/start, NULL);
+        cvReleaseImage(&img);
+        cvShowImage("avg diff", d32);
+        cvWaitKey(1);
+        printf("i: %d\r", i);
+    }
+    cvConvertScale(b32, bkg, 255, 0);
+    cvReleaseImage(&b32);
+    cvReleaseImage(&i32);
+    if (argc >= 6) {
+        cvSaveImage(argv[4], bkg, 0);
+        cvConvertScale(d32, d8, 255, 0);
+        cvSaveImage(argv[5], d8, 0); // difference image
+        return 0;
+    }
+    */
+
+    int *imgc = block_coeffs(d8, plane_coeffs);
+    IplImage *rev = splat(imgc, bsz, plane_coeffs);
+    free(imgc);
+    cvReleaseImage(&rev);
     prop_coeffs(bkg, plane_coeffs, &bkgc);
     kdt_new(&kdt, bkgc, sz, dim);
-    recon_g = cvCreateImage(bsz, bkg->depth, bkg->nChannels);
-    diff_g  = cvCreateImage(bsz, bkg->depth, bkg->nChannels);
 
-    for (i = start; i < end; i++) {
+    /*thread_ctx ctxs[3];
+    pthread_t thrs[sizeof(ctxs)/sizeof(thread_ctx)];
+    for (i = 0; i < (int)(sizeof(ctxs)/sizeof(thread_ctx)); i++) {
+        thread_ctx *ctx = &ctxs[i];
+        ctx->start = i+1;
+        ctx->end = end;
+        ctx->nb = sizeof(ctxs)/sizeof(thread_ctx);
+        ctx->path = path;
+        ctx->outfile = argc >= 7 ? argv[6] : NULL;
+        ctx->bkg = bkg;
+        ctx->diff = d8;
+        ctx->kdt = &kdt;
+        pthread_create(&thrs[i], NULL, run_thr, ctx);
+    }
+    for (i = 0; i < (int)(sizeof(ctxs)/sizeof(thread_ctx)); i++) {
+        pthread_join(thrs[i], NULL);
+    }
+    printf("all done!\n");*/
+
+    double t;
+    char outname[1024];
+    memset(outname, '\0', sizeof(outname));
+    for (i = start; i <= end; i++) {
         IplImage *img = alignedImageFrom(mkname(path, i), 8);
-        process(&kdt, bkg, img);
+        double start = get_time();
+        if (argc >= 7) snprintf(outname, sizeof(outname), "%s/bin%06d.png", argv[6], i);
+        //process(&kdt, bkg, d8, img, outname);
+        //test(img);
+        cvShowImage("image", img);
+        t += (get_time() - start);
         if ((cvWaitKey(1)&255)==27)break; // esc
         cvReleaseImage(&img);
     }
+    //free_g();
+    kdt_free(&kdt);
+    free(bkgc);
+    cvReleaseImage(&bkg);
+    cvReleaseImage(&d8);
     return 0;
 }