/**
 * CeroNoice - a shameless imitation of the ZeroNoise algorithm
 * 
 * Input: some bracketed raw files (CR2, DNG, maybe non-Canon too)
 * Output: a merged DNG file
 */
/*
 * Copyright (C) 2013 a1ex
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the
 * Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor,
 * Boston, MA  02110-1301, USA.
 */

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <limits.h>
#include <unistd.h>
#include <fcntl.h>
#include "../../src/raw.h"

#include "optmed.h"

#include "dcraw-bridge.h"
#include "exiftool-bridge.h"

/* here we only have a global raw_info */
#define save_dng(filename) save_dng(filename, &raw_info)

#define FAIL(fmt,...) { fprintf(stderr, "Error: "); fprintf(stderr, fmt, ## __VA_ARGS__); fprintf(stderr, "\n"); exit(1); }
#define CHECK(ok, fmt,...) { if (!ok) FAIL(fmt, ## __VA_ARGS__); }

#define COERCE(x,lo,hi) MAX(MIN((x),(hi)),(lo))
#define COUNT(x)        ((int)(sizeof(x)/sizeof((x)[0])))

#define MIN(a,b) \
   ({ typeof ((a)+(b)) _a = (a); \
      typeof ((a)+(b)) _b = (b); \
     _a < _b ? _a : _b; })

#define MAX(a,b) \
   ({ typeof ((a)+(b)) _a = (a); \
       typeof ((a)+(b)) _b = (b); \
     _a > _b ? _a : _b; })

#define ABS(a) \
   ({ __typeof__ (a) _a = (a); \
     _a > 0 ? _a : -_a; })

struct raw_info raw_info = {
    .api_version = 1,
    .bits_per_pixel = 16,
    .black_level = 2048,
    .white_level = 15000,
    .cfa_pattern = 0x02010100,          // Red  Green  Green  Blue
    .calibration_illuminant1 = 1,       // Daylight
};

static int black_subtract(int left_margin, int top_margin);
static int black_subtract_simple(int left_margin, int top_margin);
static int white_detect();

static inline int raw_get_pixel16(int x, int y) {
    unsigned short * buf = raw_info.buffer;
    int value = buf[x + y * raw_info.width];
    return value;
}

static inline int raw_set_pixel16(int x, int y, int value)
{
    unsigned short * buf = raw_info.buffer;
    buf[x + y * raw_info.width] = value;
    return value;
}

int raw_get_pixel(int x, int y) {
    return raw_get_pixel16(x,y);
}

/* from 14 bit to 16 bit */
int raw_get_pixel_14to16(int x, int y) {
    return (raw_get_pixel16(x,y) << 2) & 0xFFFF;
}

static int startswith(char* str, char* prefix)
{
    char* s = str;
    char* p = prefix;
    for (; *p; s++,p++)
        if (*s != *p) return 0;
    return 1;
}

static void reverse_bytes_order(char* buf, int count)
{
    unsigned short* buf16 = (unsigned short*) buf;
    int i;
    for (i = 0; i < count/2; i++)
    {
        unsigned short x = buf16[i];
        buf[2*i+1] = x;
        buf[2*i] = x >> 8;
    }
}

struct ldr
{
    unsigned short * buf;   /* original image data, unscaled, as read from the raw file */
    int black;              /* black level */
    int white;              /* white level */
    double a;               /* image matching parameters: y = ax + b */
    double b;               /* where x = in_raw - in_black, y = out_raw - out_black */
    int brightness;         /* for sorting from bright to dark */
};

struct hdr
{
#ifdef DNG_FLOAT
    float * buf;            /* binary32, scaled to 0.0 ... 1.0 */
#else
    unsigned short * buf;
#endif
} output;

struct ldr exposures[10];

static int num_exposures = 0;

static void new_exposure()
{
    if (num_exposures >= COUNT(exposures))
    {
        printf("Too many exposures (max %d)\n", COUNT(exposures));
        return;
    }

    int w = raw_info.width;
    int h = raw_info.height;
    
    unsigned short * buf = malloc(w * h * sizeof(buf[0]));
    
    int brightness = 0;
    int x,y;
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            int p = raw_get_pixel(x, y);
            brightness += p / 1024;
            buf[x + y*w] = p;
        }
    }

    exposures[num_exposures++] = (struct ldr) {
        .buf = buf,
        .black = raw_info.black_level,
        .white = raw_info.white_level,
        .a = 1,
        .b = 0,
        .brightness = brightness,
    };
}

/* call these from the darkest exposure to the brightest */
static void add_exposure(struct ldr * ldr, struct hdr * hdr)
{
    int w = raw_info.width;
    int h = raw_info.height;
    
    int x,y;
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            int raw = ldr->buf[x + y*w];
            
            /* is this pixel correctly exposed? */
            if (raw < ldr->white)
            {
                hdr->buf[x + y*w] = MAX(0, (raw - ldr->black) * ldr->a + ldr->b + raw_info.black_level);
            }
        }
    }
}

static int match_histograms(struct ldr * dark, struct ldr * bright)
{
    /* guess ISO - find the factor and the offset for matching the bright and dark images */
    int w = raw_info.width;
    int h = raw_info.height;

    int x,y;
    int i;
    
    /* build two histograms */
    int hist_size = 65536 * sizeof(int);
    int* hist_lo = malloc(hist_size);
    int* hist_hi = malloc(hist_size);
    memset(hist_lo, 0, hist_size);
    memset(hist_hi, 0, hist_size);
    
    int y0 = MAX(0, raw_info.active_area.y1 - 8) & ~3; /* gibberish above */
    
    for (y = y0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            hist_hi[bright->buf[x + y*w]]++;
            hist_lo[dark  ->buf[x + y*w]]++;
        }
    }

    /* compare the two histograms and plot the curve between the two exposures (dark as a function of bright) */
    const int min_pix = 100;                                /* extract a data point every N image pixels */
    int data_size = (w * h / min_pix + 1);                  /* max number of data points */
    int* data_x = malloc(data_size * sizeof(data_x[0]));
    int* data_y = malloc(data_size * sizeof(data_y[0]));
    double* data_w = malloc(data_size * sizeof(data_w[0]));
    int data_num = 0;
    
    int acc_lo = 0;
    int acc_hi = 0;
    int raw_lo = 0;
    int raw_hi = 0;
    int prev_acc_hi = 0;
    
    int hist_total = 0;
    for (i = 0; i < 65536; i++)
        hist_total += hist_hi[i];
    
    for (raw_hi = 0; raw_hi < bright->white; raw_hi++)
    {
        acc_hi += hist_hi[raw_hi];

        while (acc_lo < acc_hi)
        {
            acc_lo += hist_lo[raw_lo];
            raw_lo++;
        }
        
        if (raw_lo >= dark->white)
            break;
        
        if (acc_hi - prev_acc_hi > min_pix)
        {
            if (acc_hi > hist_total * 1 / 100 && acc_hi < hist_total * 99.99 / 100)    /* throw away outliers */
            {
                data_x[data_num] = raw_hi - bright->black;
                data_y[data_num] = raw_lo - dark->black;
                data_w[data_num] = (MAX(0, raw_hi - bright->black + 100));    /* points from higher brightness are cleaner */
                data_num++;
                prev_acc_hi = acc_hi;
            }
        }
    }

    /**
     * plain least squares
     * y = ax + b
     * a = (mean(xy) - mean(x)mean(y)) / (mean(x^2) - mean(x)^2)
     * b = mean(y) - a mean(x)
     */
    
    double mx = 0, my = 0, mxy = 0, mx2 = 0;
    double weight = 0;
    for (i = 0; i < data_num; i++)
    {
        mx += data_x[i] * data_w[i];
        my += data_y[i] * data_w[i];
        mxy += (double)data_x[i] * data_y[i] * data_w[i];
        mx2 += (double)data_x[i] * data_x[i] * data_w[i];
        weight += data_w[i];
    }
    mx /= weight;
    my /= weight;
    mxy /= weight;
    mx2 /= weight;
    double a = (mxy - mx*my) / (mx2 - mx*mx);
    double b = my - a * mx;

    #define BLACK_DELTA_THR 1000

    if (ABS(b) > BLACK_DELTA_THR)
    {
        /* sum ting wong */
        b = 0;
        a = (double) my / mx;
        printf("Black delta looks bad, skipping correction\n");
        goto after_black_correction;
    }

    /* apply the correction */
    bright->a = a;
    bright->b = b;

after_black_correction:
    
    if (0)
    {
        printf("Least squares   : y = %f*x + %f\n", a, b);
        FILE* f = fopen("iso-curve.m", "w");

        fprintf(f, "x = [");
        for (i = 0; i < data_num; i++)
            fprintf(f, "%d ", data_x[i]);
        fprintf(f, "];\n");
        
        fprintf(f, "y = [");
        for (i = 0; i < data_num; i++)
            fprintf(f, "%d ",data_y[i]);
        fprintf(f, "];\n");

        fprintf(f, "hl = [");
        for (i = 0; i < 65536; i++)
            fprintf(f, "%d ", hist_lo[i]);
        fprintf(f, "];\n");
        
        fprintf(f, "hh = [");
        for (i = 0; i < 65536; i++)
            fprintf(f, "%d ",hist_hi[i]);
        fprintf(f, "];\n");

        fprintf(f, "a = %f;\n", a);
        fprintf(f, "b = %f;\n", b);

        fprintf(f, "plot(x, y); hold on;\n");
        fprintf(f, "plot(x, y - b, 'g');\n");
        fprintf(f, "plot(x, a * x, 'r');\n");
        fprintf(f, "print -dpng iso-curve.png\n");
        fclose(f);
        
        if(system("octave --persist iso-curve.m"));
    }

    free(hist_lo);
    free(hist_hi);
    free(data_x);
    free(data_y);
    free(data_w);

    double factor = 1/a;
    if (factor < 1 || !isfinite(factor))
    {
        printf("Exposures not sorted?!\n");
    }
    return 1;
}

int main(int argc, char** argv)
{
    printf("\n");
    printf("CeroNoice - a clone of the ZeroNoise algorithm for merging bracketed images\n");
    printf("Usage: CeroNoice FOO.CR2 BAR.CR2 BAZ.CR2...\n");
    printf("Output: out.dng, 32-bit float\n");
    printf("\n");
    printf("  -- a1ex\n");
    printf("\n");
    
    if (argc <= 1)
        return 0;
    
    int k;
    int r;
    for (k = 1; k < argc; k++)
    {
        char* filename = argv[k];

        printf("\nInput file     : %s\n", filename);

        char dcraw_cmd[1000];
        snprintf(dcraw_cmd, sizeof(dcraw_cmd), "dcraw -v -i -t 0 \"%s\"", filename);
        FILE* t = popen(dcraw_cmd, "r");
        CHECK(t, "%s", filename);
        
        unsigned int model = get_model_id(filename);
        int exit_code = get_raw_info(model, &raw_info);

        CHECK(exit_code == 0, "RAW INFO INJECTION FAILED");

        int raw_width = 0, raw_height = 0;
        int out_width = 0, out_height = 0;
        
        char line[100];
        while (fgets(line, sizeof(line), t))
        {
            if (startswith(line, "Full size: "))
            {
                r = sscanf(line, "Full size: %d x %d\n", &raw_width, &raw_height);
                CHECK(r == 2, "sscanf");
            }
            else if (startswith(line, "Output size: "))
            {
                r = sscanf(line, "Output size: %d x %d\n", &out_width, &out_height);
                CHECK(r == 2, "sscanf");
            }
        }
        pclose(t);

        printf("Full size      : %d x %d\n", raw_width, raw_height);
        printf("Active area    : %d x %d\n", out_width, out_height);
        
        int left_margin = raw_width - out_width;
        int top_margin = raw_height - out_height;

        snprintf(dcraw_cmd, sizeof(dcraw_cmd), "dcraw -4 -E -c -t 0 \"%s\"", filename);
        FILE* fp = popen(dcraw_cmd, "r");
        CHECK(fp, "%s", filename);
        #ifdef _O_BINARY
        _setmode(_fileno(fp), _O_BINARY);
        #endif

        /* PGM read code from dcraw */
          int dim[3]={0,0,0}, comment=0, number=0, error=0, nd=0, c;

          if (fgetc(fp) != 'P' || fgetc(fp) != '5') error = 1;
          while (!error && nd < 3 && (c = fgetc(fp)) != EOF) {
            if (c == '#')  comment = 1;
            if (c == '\n') comment = 0;
            if (comment) continue;
            if (isdigit(c)) number = 1;
            if (number) {
              if (isdigit(c)) dim[nd] = dim[nd]*10 + c -'0';
              else if (isspace(c)) {
            number = 0;  nd++;
              } else error = 1;
            }
          }

        CHECK(!(error || nd < 3), "dcraw output is not a valid PGM file\n");

        int width = dim[0];
        int height = dim[1];
        CHECK(width == raw_width, "pgm width");
        CHECK(height == raw_height, "pgm height");

        void* buf = malloc(width * (height+1) * 2); /* 1 extra line for handling GBRG easier */
        int size = fread(buf, 1, width * height * 2, fp);
        CHECK(size == width * height * 2, "fread");
        pclose(fp);

        /* PGM is big endian, need to reverse it */
        reverse_bytes_order(buf, width * height * 2);

        raw_info.buffer = buf;
        
        /* did we read the PGM correctly? (right byte order etc) */
        //~ int i;
        //~ for (i = 0; i < 10; i++)
            //~ printf("%d ", raw_get_pixel16(i, 0));
        //~ printf("\n");
        
        raw_info.black_level = 2048;
        raw_info.white_level = 15000;

        raw_info.width = width;
        raw_info.height = height;
        raw_info.pitch = width * 2;
        raw_info.frame_size = raw_info.height * raw_info.pitch;

        raw_info.active_area.x1 = left_margin;
        raw_info.active_area.x2 = raw_info.width;
        raw_info.active_area.y1 = top_margin;
        raw_info.active_area.y2 = raw_info.height;
        raw_info.jpeg.x = 0;
        raw_info.jpeg.y = 0;
        raw_info.jpeg.width = raw_info.width - left_margin;
        raw_info.jpeg.height = raw_info.height - top_margin;

        white_detect();
        if (!black_subtract(left_margin, top_margin))
            printf("Black subtract didn't work\n");
        
        new_exposure();
        
        free(buf);
    }
    printf("\n");
    
    /* sort the exposures */
    int i,j;
    for (i = 0; i < num_exposures - 1; i++)
    {
        for (j = i + 1; j < num_exposures; j++)
        {
            if (exposures[i].brightness > exposures[j].brightness)
            {
                struct ldr aux = exposures[i];
                exposures[i] = exposures[j];
                exposures[j] = aux;
            }
        }
    }
    
    /* prepare for 32-bit float DNG output */
    output.buf = malloc(raw_info.width * raw_info.height * sizeof(output.buf[0]));
    raw_info.buffer = output.buf;
    
    #ifdef DNG_FLOAT
    raw_info.bits_per_pixel = 32;
    raw_info.pitch = raw_info.width * 4;
    raw_info.frame_size = raw_info.height * raw_info.pitch;
    #endif
    
    /* assume exposures are sorted from dark to bright */
    for (k = 0; k < num_exposures-1; k++)
    {
        printf("Matching images %d and %d... ", k+1, k+2); fflush(stdout);
        match_histograms(&exposures[k], &exposures[k+1]);

        double ev = log2(1/exposures[k+1].a);
        double off = exposures[k+1].b;

        /* make all scaling factors relative to first image */
        exposures[k+1].a *= exposures[k].a;
        exposures[k+1].b *= exposures[k].a;
        exposures[k+1].b += exposures[k].b;

        double ev_abs = log2(1/exposures[k+1].a);
        double off_abs = exposures[k+1].b;
        printf("%5.02fEV + %5.01f (absolute %5.02fEV + %.01f)\n", ev, off, ev_abs, off_abs);
    }

    /* cancel black offset on the brighter exposure (adjust the others instead) */
    for (k = 0; k < num_exposures; k++)
    {
        exposures[k].b -= exposures[num_exposures-1].b;
    }
    
    #ifdef DNG_FLOAT
    raw_info.black_level = 0;
    raw_info.white_level = 65536;
    #endif

    /* scale to black ... white range */
    double scale = (double)(raw_info.white_level - raw_info.black_level) / exposures[0].white;
    for (k = 0; k < num_exposures; k++)
    {
        exposures[k].a *= scale;
    }

    printf("Stacking images ");
    
    /* the first exposure is the darkest one, don't skip overexposed pixels */
    exposures[0].white = INT_MAX;
    
    for (k = 0; k < num_exposures; k++)
    {
        printf("."); fflush(stdout);
        add_exposure(&exposures[k], &output);
    }
    printf("\n");

    char out_filename[1000] = "out.dng";
    
    #ifdef DNG_FLOAT
    
    #endif

    reverse_bytes_order(raw_info.buffer, raw_info.frame_size);

    printf("Output file    : %s\n", out_filename);
    save_dng(out_filename);

    //~ copy_tags_from_source(filename, out_filename);
    return 0;
}

static int white_detect_brute_force()
{
    /* sometimes the white level is much lower than 15000; this would cause pink highlights */
    /* workaround: consider the white level as a little under the maximum pixel value from the raw file */
    /* caveat: bright and dark exposure may have different white levels, so we'll take the minimum value */
    /* side effect: if the image is not overexposed, it may get brightened a little; shouldn't hurt */
    
    /* ignore hot pixels when finding white level (at least 50 pixels should confirm it) */
    
    int white = raw_info.white_level * 5 / 6;
    int whites[8] = {white+500, white+500, white+500, white+500, white+500, white+500, white+500, white+500};
    int maxies[8] = {white, white, white, white, white, white, white, white};
    int counts[8] = {0, 0, 0, 0, 0, 0, 0, 0};

    int x,y;
    for (y = raw_info.active_area.y1; y < raw_info.active_area.y2; y ++)
    {
        for (x = raw_info.active_area.x1; x < raw_info.active_area.x2; x ++)
        {
            int pix = raw_get_pixel16(x, y);
            #define BIN_IDX ((y%4) + (x%2)*4)
            if (pix > maxies[BIN_IDX])
            {
                maxies[BIN_IDX] = pix;
                counts[BIN_IDX] = 1;
            }
            else if (pix == maxies[BIN_IDX])
            {
                counts[BIN_IDX]++;
                if (counts[BIN_IDX] > 50)
                {
                    whites[BIN_IDX] = maxies[BIN_IDX];
                }
            }
            #undef BIN_IDX
        }
    }

    //~ printf("%8d %8d %8d %8d %8d %8d %8d %8d\n", whites[0], whites[1], whites[2], whites[3], whites[4], whites[5], whites[6], whites[7]);
    //~ printf("%8d %8d %8d %8d %8d %8d %8d %8d\n", counts[0], counts[1], counts[2], counts[3], counts[4], counts[5], counts[6], counts[7]);
    
    int white1 = MIN(MIN(whites[0], whites[1]), MIN(whites[2], whites[3]));
    int white2 = MIN(MIN(whites[4], whites[5]), MIN(whites[6], whites[7]));
    white = MIN(white1, white2);
    raw_info.white_level = white - 500;
    printf("White level    : %d\n", raw_info.white_level);
    return 1;
}

static int white_detect()
{
    return white_detect_brute_force();

#if 0
    int w = raw_info.width;
    int p0 = raw_get_pixel16(0, 0);
    if (p0 < 10000) return 0;
    int x;
    for (x = 0; x < w; x++)
        if (raw_get_pixel16(x, 0) != p0)
            return 0;

    /* first line is white level, cool! */
    raw_info.white_level = p0 - 1000;       /* pink pixels at aggressive values */
    printf("White level    : %d\n", raw_info.white_level);
    return 1;
#endif
}


static int black_subtract(int left_margin, int top_margin)
{
#if 0
    reverse_bytes_order(raw_info.buffer, raw_info.frame_size);
    save_dng("untouched.dng");
#endif

    if (left_margin < 10) return 0;
    if (top_margin < 10) return 0;

    printf("Black borders  : %d left, %d top\n", left_margin, top_margin);

    int w = raw_info.width;
    int h = raw_info.height;
    
    int* vblack = malloc(h * sizeof(int));
    int* hblack = malloc(w * sizeof(int));
    int* aux = malloc(MAX(w,h) * sizeof(int));
    unsigned short * blackframe = malloc(w * h * sizeof(unsigned short));
    
    CHECK(vblack, "malloc");
    CHECK(hblack, "malloc");
    CHECK(blackframe, "malloc");

    /* data above this may be gibberish */
    int ymin = (top_margin*3/4) & ~3;

    /* estimate vertical correction for each line */
    int x,y;
    for (y = 0; y < h; y++)
    {
        int avg = 0;
        int num = 0;
        for (x = 2; x < left_margin - 8; x++)
        {
            avg += raw_get_pixel16(x, y);
            num++;
        }
        vblack[y] = avg / num;
    }
    
    /* perform some slight filtering (averaging) so we don't add noise to the image */
    for (y = 0; y < h; y++)
    {
        int y2;
        int avg = 0;
        int num = 0;
        for (y2 = y - 10*4; y2 < y + 10*4; y2 += 4)
        {
            if (y2 < ymin) continue;
            if (y2 >= h) continue;
            avg += vblack[y2];
            num++;
        }
        if (num > 0)
        {
            avg /= num;
            aux[y] = avg;
        }
        else
        {
            aux[y] = vblack[y];
        }
    }
    
    memcpy(vblack, aux, h * sizeof(vblack[0]));

    /* update the dark frame */
    for (y = 0; y < h; y++)
        for (x = 0; x < w; x++)
            blackframe[x + y*w] = vblack[y];
    
    /* estimate horizontal drift for each channel */
    int k;
    for (k = 0; k < 4; k++)
    {
        int y0 = ymin + k;
        int offset = 0;
        {
            int num = 0;
            for (y = y0; y < top_margin-4; y += 4)
            {
                offset += blackframe[y*w];
                num++;
            }
            offset /= num;
        }
        
        /* try to fix banding that repeats every 8 pixels */
        int xg;
        for (xg = 0; xg < 8; xg++)
        {
            for (x = xg; x < w; x += 8)
            {
                int num = 0;
                int avg = 0;
                for (y = y0; y < top_margin-4; y += 4)
                {
                    avg += raw_get_pixel16(x, y) - offset;
                    num++;
                }
                hblack[x] = avg / num;
            }
            
            /* perform some stronger filtering (averaging), since this data is a lot noisier */
            /* if we don't do that, we will add some strong FPN to the image */
            for (x = xg; x < w; x += 8)
            {
                int x2;
                int avg = 0;
                int num = 0;
                for (x2 = x - 1024; x2 < x + 1024; x2 += 8)
                {
                    if (x2 < 0) continue;
                    if (x2 >= w) continue;
                    avg += hblack[x2];
                    num++;
                }
                avg /= num;
                aux[x] = avg;
            }
            memcpy(hblack, aux, w * sizeof(hblack[0]));

            /* update the dark frame */
            for (y = y0; y < h; y += 4)
                for (x = xg; x < w; x += 8)
                    blackframe[x + y*w] += hblack[x];
        }
    }
    

#if 0 /* for debugging only */
    void* old_buffer = raw_info.buffer;
    raw_info.buffer = blackframe;
    reverse_bytes_order(raw_info.buffer, raw_info.frame_size);
    save_dng("black.dng");
    raw_info.buffer = old_buffer;
#endif

    /* subtract the dark frame, keeping the average black level */
    double avg_black = 0;
    for (y = top_margin; y < h; y++)
    {
        for (x = left_margin; x < w; x++)
        {
            avg_black += blackframe[x + y*w];
        }
    }
    avg_black /= (w * h);

    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            int p = raw_get_pixel16(x, y);
            int black_delta = avg_black - blackframe[x + y*w];
            p += black_delta;
            p = COERCE(p, 0, 16383);
            raw_set_pixel16(x, y, p);
        }
    }

    raw_info.black_level = round(avg_black);
    printf("Black level    : %d\n", raw_info.black_level);

#if 0
    reverse_bytes_order(raw_info.buffer, raw_info.frame_size);
    save_dng("subtracted.dng");
#endif

    free(vblack);
    free(hblack);
    free(blackframe);
    free(aux);
    return 1;
}


static int black_subtract_simple(int left_margin, int top_margin)
{
    if (left_margin < 10) return 0;
    if (top_margin < 10) return 0;
    
    int w = raw_info.width;
    int h = raw_info.height;

    /* average left bar */
    int x,y;
    long long avg = 0;
    int num = 0;
    for (y = 0; y < h; y++)
    {
        for (x = 2; x < left_margin - 8; x++)
        {
            int p = raw_get_pixel16(x, y);
            if (p > 0)
            {
                avg += p;
                num++;
            }
        }
    }
    
    int new_black = avg / num;
        
    int black_delta = raw_info.black_level - new_black;
    
    printf("Black adjust   : %d\n", (int)black_delta);
    raw_info.black_level += black_delta;
    
    return 1;
}