 /*
  * Khoros: $Id: lvspeckle.c,v 1.1 1991/05/10 15:41:54 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvspeckle.c,v 1.1 1991/05/10 15:41:54 khoros Exp $";
#endif

 /*
  * $Log: lvspeckle.c,v $
 * Revision 1.1  1991/05/10  15:41:54  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvspeckle.c
 >>>>
 >>>>      Program Name: vspeckle
 >>>>
 >>>> Date Last Updated: Tue Mar  5 22:30:07 1991 
 >>>>
 >>>>          Routines: lvspeckle - the library call for vspeckle
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvspeckle - library call for vspeckle
*
* Purpose:
*    
*    Reduce speckle noise from an image  by  using  the  Crimmins
*    algorithm.
*    
*    
* Input:
*    
*    image          VIFF image structure, must be a byte image.
*    
*    n              the number of times the user wishes  to  send  the
*                   image   through   the   speckle   removal  filter.
*                   Integer.  Limited to 20.
*    
*    
* Output:
*    
*    img2           VIFFF image structure, image with speckle  removal
*                   filter applied the proper number of times.
*    
*    
*
* Written By: Caryl Peterson   R. J. Fogler
*    
*    
****************************************************************/


/* -library_def */
int
lvspeckle(image, n)
struct xvimage *image;
int n;

/* -library_def_end */

/* -library_code */

{
   int nc,                 /* number of columns of image */
       nr,                 /*   "    "  rows    "    "   */
      image_size,          /*  nc * nr */
      buf_nc, buf_nr, buf_size,
      i,j;

   unsigned char *result; /* storage of final result */
   unsigned char *imgptr, *imgbuf;  /* the imgptr will be incremented while 
                                       the imgbuf will always point to the
                                       first pixel in the image */
   unsigned char *in_buf, *g_buf, *inptr, *g_ptr;
                           /* The image is surrounded by an edge one pixel
                              wide whose value is 0.  Thus the edges of
                              the original image are preserved.  */

   char   *program = "lvspeckle";

  /* Check type,... of image */
    if (!(propertype(program, image, VFF_TYP_1_BYTE, FALSE))) {
        (void) fprintf (stderr, "\n\n%s:   ", program);
        (void) fprintf (stderr, "lvspeckle: image must be of type byte\n");
        return (0);
    }
    if (!(proper_num_images (program, image, 1, FALSE))) {
        (void) fprintf (stderr, "\n\n%s:   ", program);
        (void) fprintf (stderr, "Can only work on files with one image\n\n");
        return (0);
    }
    if (!(proper_num_bands (program, image, 1, FALSE))) {
        (void) fprintf (stderr,"\n\n%s:   ", program);
        (void) fprintf (stderr,"Can only work on images with 1 data band\n\n");
        return (0);
    }


    nc = image->row_size; /* number of columns */
    nr = image->col_size; /* number of rows    */

    imgbuf = (unsigned char *)(image->imagedata); /* pointer to image's pixels */


/* *********************************************************** *
 * *        A L L O C A T E   I M A G E   B U F F E R S      * *
 * *********************************************************** */
    buf_nc = nc + 2;
    buf_nr = nr + 2;
    buf_size = buf_nc * buf_nr;

    /*  create in_buffer */
    in_buf = (unsigned char *) malloc (buf_size);
    if (in_buf == NULL) {
      (void) fprintf (stderr, "\nvspeckle unable to allocate buffer\n");
      exit (1);
     }

   /* create g buffer */
    g_buf = (unsigned char *) malloc (buf_size);
    if (g_buf == NULL) {
      (void) fprintf (stderr, "\nvspeckle unable to allocate buffer\n");
      exit (1);
     }

/* *********************************************************** *
 * *     P U T   S O U R C E   I M A G E  I N  B U F F E R   * *
 * *********************************************************** */
  imgptr = imgbuf;
  inptr = in_buf;
         /*  zero first row of buffer */
  for (i = 0; i < buf_nc; i++) {
      *inptr++ = ( unsigned char ) 0;
  }
  
  for (i = 0; i < nr ; i++) {
    *inptr++ = (unsigned char) 0;   /* zero first column */
    for (j = 0; j < nc; j++) {
       *inptr++ = *imgptr++;  /* place the image pixel's values in the buffer */
    }
    *inptr++ = (unsigned char) 0;   /* zero last column*/
     
    }
  
  for (i = 0; i < buf_nc; i++) { 
      *inptr++ = ( unsigned char ) 0;    /*  zero last row */
  }

/* *********************************************************** *
 *     P R E P A R E  G _ B U F F E R                     * *
 * *********************************************************** */
  g_ptr = g_buf;
         /*  zero first row of buffer */
  for (i = 0; i < buf_nc; i++) {
      *g_ptr++ = ( unsigned char ) 0;
  }
  for (i = 0; i < nr ; i++) {
    *g_ptr++ = (unsigned char) 0;   /* zero first column */
     g_ptr = (unsigned char *)(g_ptr + nc);
    *g_ptr++ = (unsigned char) 0;   /* zero last column*/
    }

  for (i = 0; i < buf_nc; i++) {  
      *g_ptr++ = ( unsigned char ) 0;  /* zero last row */
  }


/* *********************************************************** *
 * *    S P E C K L E   R E M O V A L   A L G O R I T H M    * *
 * *********************************************************** */


  for (j = 0; j < n; j++) {
     geomfilt(nr, nc, in_buf, g_buf);
   }

/* *********************************************************** *
 * *    P R E P A R E   R E S U L T A N T    I M A G E       * *
 * *********************************************************** */

 imgptr = imgbuf;
 inptr = in_buf;
 
 inptr = (unsigned char *) ( inptr + buf_nc);  /* skip first row*/

 for (i = 0; i < nr; i++){

   inptr++;  /* skip first column */
   
    for (j = 0; j < nc; j++) {
      *imgptr++ = *inptr++;
    }
   inptr++;   /* skip last column */
  }
  return (1);
}

/* *********************************************************** *
 * *   G E O M E T R I C   F I L T E R   P R O C E D U R E   * *
 * *********************************************************** */
int geomfilt(M, N, f, g)
int M, N;
unsigned char *f, *g;
{
  symhull(1, 0, M, N, f, g);   /* N-S   symmetrical hull */
  symhull(0, 1, M, N, f, g);   /* E-W   symmetrical hull */
  symhull(1, 1, M, N, f, g);   /* NE-SW symmetrical hull */
  symhull(1,-1, M, N, f, g);   /* NW-SE symmetrical hull */
}



/* *********************************************************** *
 * *    S Y M M E T R I C A L   H U L L   P R O C E D U R E  * *
 * *********************************************************** */
int symhull(a, b, M, N, f, g)
int a, b, M, N;
unsigned char *f, *g;
{
          /* send the image twice to the positive hull */
  poshull( a, b, M, N, f, g);
  poshull(-a,-b, M, N, f, g);
          /* send the image twice to the negative hull */
  neghull(-a,-b, M, N, f, g);
  neghull( a, b, M, N, f, g);

}


/* *********************************************************** *
 * *       P O S I T I V E   H U L L   P R O C E D U R E     * *
 * *********************************************************** */
int poshull(a, b, M, N, f, g)
int a, b, M, N;
unsigned char *f, *g;
{
  int m, n;
  int mn, ab, ba;
  unsigned char tmp, temp1, temp2;
  unsigned char *fmn, *gmn, *fab, *gab, *gba;
  unsigned char *fmn_p, *gmn_p, *fab_p, *gab_p, *gba_p;

  mn = (N+2) + 1;
  ba = a*(N+2) + b;
  ab = -ba;
  fmn = (unsigned char *)(f + mn);
  gmn = (unsigned char *)(g + mn);
  fab = (unsigned char *)(fmn + ab);
  gab = (unsigned char *)(gmn + ab);
  gba = (unsigned char *)(gmn + ba);

  fmn_p = fmn;
  fab_p = fab;
  gmn_p = gmn;
  for (m = 1 ; m <= M ; m++) {
    for (n = 1 ; n <= N ; n++) {
      tmp = *fmn_p++;   /* this is the middle pixel */
      if ((*fab_p) < (tmp +1)) {
         temp1 = (*fab_p);
       }
      else{
         temp1 = tmp+1;
       }
       fab_p++;
       if (tmp > temp1) {
           *gmn_p++ = tmp;
        }
        else {
           *gmn_p++ = temp1;
         }
    }
    fmn_p++;
    gmn_p++;
    fab_p++;
  }

  fmn_p = fmn;
  gmn_p = gmn;
  gab_p = gab;
  gba_p = gba;
  for (m = 1 ; m <= M ; m++) {
    for (n = 1 ; n <= N ; n++) {
      tmp = *gmn_p++;  /* the middle pixel */
      if ((tmp + 1) < (*gba_p + 1 ) ) {
          temp1 = tmp + 1; 
        }
      else{
          temp1 = *gba_p + 1;
        }
      gba_p++;
      if ((*gab_p) < temp1) {
           temp2 = *gab_p;
       }
      else {
           temp2 = temp1;
       }
      gab_p++;
      if (tmp > temp2 ) {
          *fmn_p++ = tmp;
       }
      else {
          *fmn_p++ = temp2;
       }
    }
    fmn_p++;
    gmn_p++;
    gab_p++;
    gba_p++;
  }
}

/* *********************************************************** *
 * *      N E G A T I V E   H U L L   P R O C E D U R E      * *
 * *********************************************************** */
int neghull(a, b, M, N, f, g)
int a, b, M, N;
unsigned char *f, *g;
{
  int m, n;
  int mn, ab, ba;
  unsigned char tmp,  temp1, temp2;
  unsigned char *fmn, *gmn, *fab, *gab, *gba;
  unsigned char *fmn_p, *gmn_p, *fab_p, *gab_p, *gba_p;

  mn = (N+2) + 1;
  ba = a*(N+2) + b;
  ab = -ba;
  fmn = (unsigned char *)(f + mn);
  gmn = (unsigned char *)(g + mn);
  fab = (unsigned char *)(fmn + ab);
  gab = (unsigned char *)(gmn + ab);
  gba = (unsigned char *)(gmn + ba);

  fmn_p = fmn;
  gmn_p = gmn;
  fab_p = fab;
  for (m = 1 ; m <= M ; m++) {
    for (n = 1 ; n <= N ; n++) {
      tmp = *fmn_p++;  /* the middle pixel */
      if( ((*fab_p) +1) > (tmp - 1)) {
          temp1 = (*fab_p) + 1;
       }
      else {
          temp1 = tmp - 1;
      }
      fab_p++;
      if (tmp < temp1){
         *gmn_p++ = tmp;
       }
      else {
          *gmn_p++ = temp1;
       }
       
    }
    fmn_p++;
    gmn_p++;
    fab_p++;
  }

  fmn_p = fmn;
  gmn_p = gmn;
  gab_p = gab;
  gba_p = gba;
  for (m = 1 ; m <= M ; m++) {
    for (n = 1 ; n <= N ; n++) {
      tmp = *gmn_p++;  /* middle pixel */
      if ((tmp-1) > ((*gba_p) +1)) {
         temp1 = tmp-1;
       }
      else {
         temp1 = (*gba_p) +1;
       }
      gba_p++;

      if (*gab_p > temp1) {
          temp2 = *gab_p;
       }
      else {
           temp2 = temp1;
       }
      gab_p++;
      if (tmp < temp2) {
          *fmn_p++ = tmp;
       }
      else {
           *fmn_p++ = temp2;
       }
    }

    fmn_p++;
    gmn_p++;
    gab_p++;
    gba_p++;
  }
}


/* -library_code_end */
