/* See {dsm_segment.h}. */
/* Last edited on 2016-04-09 11:03:55 by stolfilocal */

#define _GNU_SOURCE
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

#include <dsm_image.h>
#include <dsm_segment.h>

void dsm_segment_compute_mismatch
    ( uint8_t *img, 
      uint8_t *img_blur,
      uint8_t *bg,
      uint8_t *bg_blur,
      int nx,
      int ny, 
      int x, 
      int y,
      double *d, double *w )
  {
    /*Compute the local average of img and bg.*/
    double avg_p = (double)img_blur[y * nx + x];
    double avg_b = (double)bg_blur[y * nx + x];
    double ctr_p = (double)img[y * nx + x];
    double ctr_b = (double)bg[y * nx + x];
    double avg_coeff = (avg_p + 1) / (avg_b + 1);
    double ctr_coeff = (ctr_p + 1) / (ctr_b + 1);
    *w = fmin (avg_p, avg_b);
    *d = log(ctr_coeff/avg_coeff)/M_LN2;
  }

  void dsm_segment_identify_background
    ( dsm_image_t *img, 
      dsm_image_t *img_blur,
      dsm_image_t *bg,
      dsm_image_t *bg_blur, 
      double lambda,
      dsm_image_t *classw
    )
  {
    int nx = img->nx;
    int ny = img->ny;
    assert((nx == img_blur->nx) && (ny == img_blur->ny));
    assert((nx == bg->nx) && (ny == bg->ny));
    assert((nx == bg_blur->nx) && (ny == bg_blur->ny));
    for (int y = 0; y < ny; y++) {
      for (int x = 0; x < nx; x++) {
	double d, w;
	dsm_segment_compute_mismatch (img->p, img_blur->p, bg->p, bg_blur->p, nx, ny, x, y, &d, &w);
	double z = d/lambda;
	double cw = exp (-z*z/2);
	classw->p[y * nx + x] = dsm_image_quantize(cw * 255, 0.0, 255.0); 
      }
    }
  }