 /*
  * Khoros: $Id: lvgettilt.c,v 1.2 1991/10/02 00:03:37 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvgettilt.c,v 1.2 1991/10/02 00:03:37 khoros Exp $";
#endif

 /*
  * $Log: lvgettilt.c,v $
 * Revision 1.2  1991/10/02  00:03:37  khoros
 * HellPatch2
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * 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 to 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: lvgettilt.c
 >>>>
 >>>>      Program Name: vgettilt
 >>>>
 >>>> Date Last Updated: Thu Sep 26 11:12:07 1991 
 >>>>
 >>>>          Routines: lvgettilt - the library call for vgettilt
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: lvgettilt - library call for vgettilt
*
* Purpose:
*    
*    computes the best-fit plane parameters for  the  input  VIFF
*    file
*    
*    

* Input:
*    
*    lvgettilt:
*    
*    image          (struct xvimage) - input xvimage structure
*    
*    mimage         (struct xvimage) - input mask xvimage structure
*    
*    mflg           (int) - flag set to 1 if masking image  is  speci-
*                   fied
*    
*                   lvfindtilt:
*    
*    image          (struct xvimage) - input xvimage structure
*    
*                   lvmafdtilt:
*    
*    image          (struct xvimage) - input xvimage structure
*    
*    mimage         (struct xvimage) - input mask xvimage structure
*    
*    

* Output:
*    
*    lvgettilt, lvfindtilt, & lvmafdtilt:
*    
*    a              (double) - A coefficient of f(x,y)=Ax+By+C
*    
*    b              (double) - B coefficient of f(x,y)=Ax+By+C
*    
*    c              (double) - C coefficient of f(x,y)=Ax+By+C
*    
*    These library functions will return a zero if failure and  a  one
*    upon success.
*    
*    

*
* Written By: Scott Wilson, Donna Koechner
*    
*    Mon Jun 11 16:30:29 MDT 1990 Donna Koechner - updated  to  ghost-
*    writer
*    
*    Wed Mar 20 12:36:43 MET 1991 Jean Pierre Guerin
*    
*             1- fixed loops to work right with non square images
*    
*             2- fixed bogus (int) caste
*    
*             3- fixed incorrect term in calculation of detc
*    
*             4- corrected pointer increments
*    
*    25-Sep-91 Scott Wilson - Changed (int) castes to  (unsigned  int)
*    castes  in  byte  processing  loop - otherwise can get goofy sign
*    extension problems.
*    
*    25-Sep-91 Scott Wilson - Found that the det{A,B,C} eqns were  ALL
*    wrong.  Fixed them up after rederiving the whole mess.
*    
*    26-Sep-91 Scott Wilson - Fixed incorrect  center  point  calcula-
*    tion.
*    
*    

****************************************************************/


/* -library_def */
lvgettilt(image, mimage, mflg, a, b, c)
struct xvimage  *image, *mimage;
int             mflg;
double          *a,*b,*c;
/* -library_def_end */

/* -library_code */
{
  if (mflg)
  {
    if (!lvmafdtilt(image, mimage, a, b, c) )
    {
      (void) fprintf(stderr, "lvmafdtilt Failed\n");
      return(FALSE);
    }
  }
  else  /* no masking image */
  {
    if (! lvfindtilt(image, a, b, c))
    {
      (void) fprintf(stderr, "lvfindtilt Failed\n");
      return(FALSE);
    }
  }

  return(TRUE);

}

/**************************************************************
*
* MODULE NAME: lvfindtilt
*
*     PURPOSE: Compute the least-squares best fit plane for an image.
#                The plane will be given by f(x,y)=Ax+By+C
*                with  the tilt coefficients in inverse pixels.
*
*       INPUT:
*             image - input xvimage structure
*
*      OUTPUT: 
*              a,b,c - specify the coefficients of the plane.
*
*              lvfindtilt will return a zero if failure and a
*              one upon success.
*
* CALLED FROM: main
*
* ROUTINES CALLED:
*
**************************************************************/

int
lvfindtilt(image,a,b,c)
struct xvimage *image;
double *a,*b,*c;
  {
    double x2,xy,x,y2,y,xz,yz,z,n,det,deta,detb,detc, i_rxi_r, i_r,j_c;
    double rctr,cctr;
    int rows,cols;
    int i,j;
    unsigned char *ch;
    short *s;
    int *it;
    float *f;

    x2 = 0.0; xy = 0.0; x = 0.0; y2 = 0.0; y = 0.0; xz = 0.0;
    yz = 0.0; z  = 0.0;

    rows = image->col_size;
    cols = image->row_size;
    rctr = (rows-1)/2.0;
    cctr = (cols-1)/2.0;
    if (rctr == 0 || cctr == 0)
      {
        (void) fprintf(stderr,"lvfindtilt: Cannot find tilt on one pixel wide image!\n");
        return(0);
      }

    n = rows*cols;

    switch(image->data_storage_type)
      {
        case VFF_TYP_BIT:
          (void) fprintf(stderr,"lvfindtilt: Cannot find tilt for BIT image!\n");
          return(0);
          break;
        case VFF_TYP_1_BYTE:
          ch = (unsigned char *)image->imagedata;    
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, ch++)
                {
                  j_c= j - cctr;
                  x  += j_c;
                  x2 += j_c * j_c;
                  xy += j_c * i_r;
                  y  += i_r;
                  y2 += i_rxi_r;
                  xz += j_c * (unsigned int)*ch;
                  yz += i_r * (unsigned int)*ch;
                  z  += (unsigned int)*ch;
                }
            }
          break;
        case VFF_TYP_2_BYTE:
          s = (short *)image->imagedata;    
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, s++)
                {
                  j_c= j - cctr;
                  x  += j_c;
                  x2 += j_c * j_c;
                  xy += j_c * i_r;
                  y  += i_r;
                  y2 += i_rxi_r;
                  xz += j_c * (*s);
                  yz += i_r * (*s);
                  z  += *s;
                }
            }
          break;
        case VFF_TYP_4_BYTE:
          it = (int *)image->imagedata;    
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, it++)
                {
                  j_c= j - cctr;
                  x  += j_c;
                  x2 += j_c * j_c;
                  xy += j_c * i_r;
                  y  += i_r;
                  y2 += i_rxi_r;
                  xz += j_c * (*it);
                  yz += i_r * (*it);
                  z  += *it;
                }
            }
          break;
        case VFF_TYP_FLOAT:
          f = (float *)image->imagedata;    
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, f++)
                {
                  j_c= j - cctr;
                  x  += j_c;
                  x2 += j_c * j_c;
                  xy += j_c * i_r;
                  y  += i_r;
                  y2 += i_rxi_r;
                  xz += j_c * (*f);
                  yz += i_r * (*f);
                  z  += *f;
                }
            }
          break;
        case VFF_TYP_COMPLEX:
          (void) fprintf(stderr,"lvfindtilt: Cannot find tilt for COMPLEX image!\n");
          return(0);
          break;
        default:
          (void) fprintf(stderr,"lvfindtilt: Unkown data storage type: %d\n",image->data_storage_type);
          return(0);
          break;
      }

    /* Find the matrix determinant */
    det = x2*(y2*n-y*y)-xy*(xy*n-y*x)+x*(xy*y-y2*x);

    /* Find the Cramer's determinants */
    deta = xz*(y2*n-y*y)-xy*(yz*n-y*z)+x*(yz*y-y2*z);
    detb = x2*(yz*n-y*z)-xz*(xy*n-x*y)+x*(xy*z-yz*x);
    detc = x2*(y2*z-yz*y)-xy*(xy*z-yz*x)+xz*(xy*y-y2*x);

    /* Compute the plane fit coefficients */
    *a = deta/det;
    *b = detb/det;
    *c = detc/det;

    return(1);  /* success */

  }


/**************************************************************
*
* MODULE NAME: lvmafdtilt
*
*     PURPOSE: Compute the least-squares best fit plane for an image.
*                The plane will be given by f(x,y)=Ax+By+Cz
*                with  the tilt coefficients in inverse pixels.
*                The plane coefficients will be computed only for those
*                pixels enabled by the mask image.
*
*       INPUT:
*             image - input xvimage structure
*             mimage - input mask xvimage structure
*
*             The mask image and the input image MUST be of the same
*             type and size!
*
*      OUTPUT: 
*              a,b,c - specify the coefficients of the plane.
*
*              lvmafdtilt will return a zero if failure and a
*              one upon success.
*
* CALLED FROM: main
*
* ROUTINES CALLED:
*
**************************************************************/

int
lvmafdtilt(image,mimage,a,b,c)
struct xvimage *image,*mimage;
double *a,*b,*c;
  {
    double x2,xy,x,y2,y,xz,yz,z,n,det,deta,detb,detc, i_rxi_r, i_r,j_c;
    double rctr,cctr;
    int rows,cols;
    int i,j;
    unsigned char *ch,*mch;
    short *s,*ms;
    int *it,*mit;
    float *f,*mf;

    x2 = 0.0; xy = 0.0; x = 0.0; y2 = 0.0; y = 0.0; xz = 0.0;
    yz = 0.0; z  = 0.0;

    if (image->data_storage_type != mimage->data_storage_type)
      {
        (void) fprintf(stderr,"lvmafdtilt: Mismatched input and mask!\n");
        return(0);
      }
    rows = image->col_size;
    cols = image->row_size;
    rctr = (rows-1)/2.0;
    cctr = (cols-1)/2.0;
    if (rctr == 0 || cctr == 0)
      {
        (void) fprintf(stderr,"lvmafdtilt: Cannot find tilt for a one pixel wide image!\n");
        return(0);
      }

    n = 0;

    switch(image->data_storage_type)
      {
        case VFF_TYP_BIT:
          (void) fprintf(stderr,"lvmafdtilt: Cannot find tilt for BIT image!\n");
          return(0);
          break;
        case VFF_TYP_1_BYTE:
          ch = (unsigned char *)image->imagedata;    
          mch = (unsigned char *)mimage->imagedata;
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, ch++)
                {
                  if (*mch++)
                    {
                      j_c= j - cctr;
                      x  += j_c;
                      x2 += j_c * j_c;
                      xy += j_c * i_r;
                      y  += i_r;
                      y2 += i_rxi_r;
                      xz += j_c * (unsigned int)*ch;
                      yz += i_r * (unsigned int)*ch;
                      z  += (unsigned int)*ch;
                      n++;
                    }
                }
            }
          break;
        case VFF_TYP_2_BYTE:
          s = (short *)image->imagedata;    
          ms = (short *)mimage->imagedata;
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, s++)
                {
                  if (*ms++)
                    {
                      j_c= j - cctr;
                      x  += j_c;
                      x2 += j_c * j_c;
                      xy += j_c * i_r;
                      y  += i_r;
                      y2 += i_rxi_r;
                      xz += j_c * (*s);
                      yz += i_r * (*s);
                      z  += *s;
                      n++;
                    }
                }
            }
          break;
        case VFF_TYP_4_BYTE:
          it = (int *)image->imagedata;    
          mit = (int *)mimage->imagedata;
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, it++, mit++)
                {                                 /* i++ replaced by it++ */
                  if (*mit)
                    {
                      j_c = j - cctr;
                      x  += j_c;
                      x2 += j_c * j_c;
                      xy += j_c * i_r;
                      y  += i_r;
                      y2 += i_rxi_r;
                      xz += j_c * (*it);
                      yz += i_r * (*it);
                      z  += (*it);
                      n++;
                    }
                }
            }
          break;
        case VFF_TYP_FLOAT:
          f = (float *)image->imagedata;    
          mf = (float *)mimage->imagedata;
          for (i=0; i<rows; i++)
            { 
              i_r     = i - rctr;
              i_rxi_r = i_r * i_r;
              for (j=0; j<cols; j++, f++)
                {
                  if (*mf++)
                    {
                      j_c= j - cctr;
                      x  += j_c;
                      x2 += j_c * j_c;
                      xy += j_c * i_r;
                      y  += i_r;
                      y2 += i_rxi_r;
                      xz += j_c * (*f);
                      yz += i_r * (*f);
                      z  += *f;
                      n++;
                    }
                }
            }
          break;
        case VFF_TYP_COMPLEX:
          (void) fprintf(stderr,"lvmafdtilt: Cannot find tilt for COMPELX image!\n");
          return(0);
          break;
        default:
          (void) fprintf(stderr,"lvmafdtilt: Unknown data storage type: %d\n",image->data_storage_type);
          return(0);
          break;
      }

    /* Find the matrix determinant */
    det = x2*(y2*n-y*y)-xy*(xy*n-y*x)+x*(xy*y-y2*x);

    /* Find the Cramer's determinants */
    deta = xz*(y2*n-y*y)-xy*(yz*n-y*z)+x*(yz*y-y2*z);
    detb = x2*(yz*n-y*z)-xz*(xy*n-x*y)+x*(xy*z-yz*x);
    detc = x2*(y2*z-yz*y)-xy*(xy*z-yz*x)+xz*(xy*y-y2*x);

    /* Compute the plane fit coefficients */
    *a = deta/det;
    *b = detb/det;
    *c = detc/det;

    return(1);  /* success */

  }

/* -library_code_end */
