/*Dezing:	This routine dezingers a pair of images, putting the output
                in the first image.

Derived from the ADSC routine dezinger_to_int_simple, March 2003
*/


#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "objects.h"

void dezing(unsigned short *indat1,unsigned short *indat2, int slow_max,
 int fast_max, int saturated, float *tempbuf)
{
  unsigned short *ptr1, *ptr2;
  float *fptr;
  int i, j, l, ifm, ifn, n, diff_sig, pix1, pix2, npix, linct, mini;
  float avedif, rhold;

  diff_sig = 8;
  npix = slow_max * fast_max;
  linct = 100 * fast_max;
  printf ("Dezingering:");

  ptr1 = indat1;
  ptr2 = indat2;
  ifn = 0;
  mini = 65536;
  n = 0;

  for (i = 0; i < npix; i++) {
    pix1 = *ptr1;
    pix2 = *ptr2;
    if ((pix1 < saturated) && (pix2 < saturated)) {
      ifn += abs(pix1 - pix2);
      n++;
    }
    if ((pix1 < mini) && (pix1 > 0))
      mini = pix1;
    ptr1++;
    ptr2++;
  }
  avedif = ifn;
  rhold = 2*n;
  if (n > 0)
    avedif = avedif/rhold;
  else
    avedif = 0.0;
  printf (" avg abs diff (%d pix) %6.2f,",n,avedif);

  l = avedif;
  l = diff_sig * (1 + l);

  ifn = 0;
  j = 0;
  ptr1 = indat1;
  ptr2 = indat2;
  for (i = 0; i < npix; i++) {
    pix1 = *ptr1;
    pix2 = *ptr2;
    if ((pix1 < saturated) && (pix2 < saturated)) {
      ifm = abs(pix1 - pix2);
      if (ifm < l) {
        ifn += ifm;
        j++;
      }
    }
    ptr1++;
    ptr2++;
  }
  avedif = ifn;
  rhold = 2*j;
  if (j > 0)
    avedif = avedif/rhold;
  else
    avedif = 0.0;
  printf (" after rejs %6.2f\n",avedif);

  l = avedif;
  l = diff_sig * (1 + l);

  n = 0;
  ptr1 = indat1;
  ptr2 = indat2;
  fptr = tempbuf;
  for (i = 0; i < npix; i++) {
    pix1 = *ptr1;
    pix2 = *ptr2;
    ifm = abs(pix1 - pix2);
    if ((ifm < l) && (pix1 < saturated) && (pix2 < saturated)) {
      rhold = pix1 + pix2;            /* Both pixels unsat'd and well-behaved */
      *fptr = 0.5 * rhold;
    }
    else {          /* difference is too large, one is a zinger or saturated */
      n++;
      if (pix2 < saturated) {               /* 2nd exposure un-saturated */
        if (pix1 < pix2) {                  /* 1st exposure also */
          *fptr = pix1;                     /* Use first */
        }
        else {                              /* Second smaller, use it */
          *fptr = pix2;
        }
      }
      else {                               /* 2nd point saturated, use first */
        *fptr = pix1;
      }
    }
    ptr1++;
    ptr2++;
    fptr++;
    if (i%linct == 0) {
      j = i/fast_max;
      printf("%d ",j);
      fflush(stdout);
    }
  }
  printf (" - dezingered %d pixels\n",n);
  ptr1 = indat1;
  fptr = tempbuf;
  for (i = 0; i < npix; i++) {    /* Copy data back to original array */
    pix1 = *fptr;
    if (pix1 > 65535)
      pix1 = 65535;
    if (pix1 < 0)
      pix1 = 0;
    *ptr1 = pix1;
    ptr1++;
    fptr++;
  }

  printf("\n");
  return;
}