 /*
  * Khoros: $Id: ldmpp.c,v 1.3 1992/03/20 23:28:27 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldmpp.c,v 1.3 1992/03/20 23:28:27 dkhoros Exp $";
#endif

 /*
  * $Log: ldmpp.c,v $
 * Revision 1.3  1992/03/20  23:28:27  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, 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 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldmpp.c
 >>>>
 >>>>      Program Name: dmpp
 >>>>
 >>>> Date Last Updated: Mon Mar  9 18:54:00 1992 
 >>>>
 >>>>          Routines: ldmpp - the library call for dmpp
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define CONSTANT 57.295779
static int ldmag(),ldpower(),ldphase();
/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldmpp - library call for dmpp
*
* Purpose:
*    
*    Magnitude, Power, or Phase of 1-D signal(s)
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    mode           determines type of data desired (0 = magnitude,  1
*                   = power, 2 = phase).
*    
*    output         magnitude or power output option:  0  =  magnitude
*                   or  power;  1  =  log(magnitude  or  power);  2  =
*                   log(magnitude or power +1); 3  =  ln(magnitude  or
*                   power); 4 = ln(magnitude or power +1).
*    
*    phase          phase output option:  0 = radians; 1 = degrees.
*    
*    direction      process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

* Output:
*    
*    image          pointer to VIFF structure  containing  image  data
*                   after processing.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Ramiro Jordan, Jeremy Worley
*    
*    Jeremy Worley 19 Dec 1990 18:18 MST
*              Modified phase calculation to explicitly check for  the
*              (0,0)  and  set  the  phase at those points to 0.  This
*              prevents the dreaded "SQRT domain error"
*    
*    Jeremy Worley 05 Mar 1992 10:52 MST
*              Added a missing declaration for the argument output  in
*              ldmpp().  Added an explicit declaration of propertype()
*              to ldmag().
*    
*    

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


/* -library_def */
int
ldmpp (image,mode,output,phase,min_phase,direction)
   struct xvimage  *image;
   int             output,mode,phase,min_phase,direction;
/* -library_def_end */

/* -library_code */
{
  char *program = "ldmpp";

  switch(mode){
    case 0 :  if(!ldmag(image,output,direction)){
                 fprintf(stderr,"%s:  ldmag failed.\n",program);
                 return(0);
              }
              break;   
    case 1 :  if(!ldpower(image,output,direction)){
                 fprintf(stderr,"%s:  ldpower failed.\n",program);
                 return(0);
              }
              break;
    case 2 :  if(!ldphase(image,phase,min_phase,direction)){
                 fprintf(stderr,"%s:  ldphase failed.\n",program);
                 return(0);
              }
              break;
    default : fprintf(stderr,"%s:  illegal mode specified.\n",program);
              return(0);
  }

  return(1);
} 

/***********************************************************************
*
*  Routine Name: ldmag
*
*          Date: May 24, 1990
*        
*       Purpose: Calculates magnitude of 1-D signal(s) 
*
*         Input: *image         pointer to the input/output xvimage 
*                               structure
*
*                magnitude_type 0 = magnitude, 1 = log(magnitude),
*                               2 = log(magnitude+1), 3 = ln(magnitude),
*                               4 = ln(magnitude+1).
*
*                process_dir    (int) data processing direction.  
*                               0 = vector, 1 = band.
*
*        Output: *image         (struct xvimage) - pointer to the input/output
*                               xvimage structure.
*
*    Written By: Ramiro Jordan 
*
* Modifications: Jeremy Worley 14-Aug-90 Incorporated into dmpp and made
*                                        static.
*
***********************************************************************/

static int ldmag ( image, magnitude_type, process_dir)
struct xvimage  *image;
int             magnitude_type,  /* magnitude type */
                process_dir;     /* direction of processing 0=vect, 1= band */
{
    int    i, j, k, points;
    int    dunload_vector();
    int    dimension, numvects;
    float  **vector;
    char   **dload_vector();
    double cmagsq();
    int propertype();

    char *program = "ldmpp";

/*
** perform error checking 
*/

    if(!(propertype(program, image, VFF_TYP_FLOAT, FALSE)) &&
       !(propertype(program, image, VFF_TYP_COMPLEX, FALSE))) {
       fprintf(stderr,"\n\n%s:  ", program);
       fprintf(stderr,"data storage type must be of type float or complex\n");
       return (0);
    }

/*
** Call dload_vector for input image
*/

    if((vector=(float **)dload_vector(image,&numvects,&dimension,process_dir))
        == NULL) {
       fprintf(stderr,"%s:  Unable to load vector.\n", program);
       return(0);
    }

/*
** Perform the desired operation on the data sequence
*/

   for(k=0;k<numvects;k++) {
       if(image->data_storage_type == VFF_TYP_FLOAT) {

          for(i=0;i<dimension;i++) {
              switch (magnitude_type) {
                 case 0:  break;
                 case 1:  vector[k][i] = log10((double)vector[k][i]);
                          break;
                 case 2:  vector[k][i] = log10((double)vector[k][i] + 1.0);
                          break;
                 case 3:  vector[k][i] = log((double)vector[k][i]);
                          break;
                 case 4:  vector[k][i] = log((double)vector[k][i] + 1.0);
                          break;
                 default: break;
              } /* end switch */
          } /* end for */

       } else if(image->data_storage_type == VFF_TYP_COMPLEX){
          points = dimension * 2;

          for(i=0,j=0;i<points;i+=2,j++){
              vector[k][j] = sqrt((double)cmagsq(vector[k][i],vector[k][i+1]));
              switch(magnitude_type){
                 case 0:  break;
                 case 1:  vector[k][j] = log10((double)vector[k][j]);
                          break;
                 case 2:  vector[k][j] = log10((double)vector[k][j] + 1.0);
                          break;
                 case 3:  vector[k][j] = log((double)vector[k][j]);
                          break;
                 case 4:  vector[k][j] = log((double)vector[k][j] + 1.0);
                          break;
                 default: break;
              } /* end switch */
          } /* end for */
       } /* end big if */
    } /* end for */

/*
** Call dunload to return vector to its "original" format
*/

    if(!dunload_vector((char **)vector,image,(unsigned long)VFF_TYP_FLOAT,
        numvects,dimension, process_dir)) {
       fprintf(stderr,"%s:  Unable to unload vector.\n", program);
       return(0);
    }

    for(i=0;i<numvects;i++)free(vector[i]);
    free(vector); 

    return(1);
} /* end of lmag */

/***********************************************************************
*
*  Routine Name: ldpower()
*
*          Date: July 12, 1990 (?)
*        
*       Purpose: Power of 1-D signal(s) 
*
*         Input: image          pointer to VIFF structure containing image 
*                               data to be processed.
*
*                power_type     power output option:  0 = magnitude, 
*                               1 = log(power), 2 = log(power+1), 
*                               3 = ln(power), 4 = ln(power+1).
*
*                procdir        process direction:  0 indicated vector 
*                               oriented processing, 1 indicates band 
*                               oriented processing.
*
*        Output: image          pointer to VIFF structure containing image
*                               data after processing.
*
*    Written By: Ramiro Jordan 
*
* Modifications: Jeremy Worley 14-Aug-90 Incorporated into dmpp and made
*                                        static.
*
***********************************************************************/

static int ldpower ( image, power_type, process_dir)

struct xvimage  *image;
int             power_type,      /* power type */
                process_dir;     /* direction of processing 0=vect, 1= band */
{
    int    i, j, k, points;
    int    dunload_vector();
    int    dimension, numvects;
    float  **vector;
    char   **dload_vector();
    double cmagsq();

    char *program = "ldpower";

    /* perform error checking */
    if ( !(propertype(program, image, VFF_TYP_FLOAT, FALSE)) &&
         !(propertype(program, image, VFF_TYP_COMPLEX, FALSE)) )
    {
          (void) fprintf (stderr, "\n\n%s:   ", program);
          (void) fprintf (stderr, "ldpower: data storage type must be of type, float or complex\n\n");
          return (0);
    }

/*
 * Call dload_vector for input image
 */
    if ( (vector = (float **)dload_vector(image,&numvects,&dimension,process_dir)) == NULL )
    {
       fprintf(stderr,"%s:  Unable to load vector.\n", program);
       return(0);
    }

/*
 * Perform the desired operation on the data sequence
 */
   for ( k = 0; k < numvects; k++ )
   {
        if ( image->data_storage_type == VFF_TYP_FLOAT )
        {
           for ( i = 0; i < dimension; i++ )
           {
               vector[k][i] *=  vector[k][i];
               switch ( power_type )
               {
                   case 0:
                      break;
                   case 1:
                      vector[k][i] = (float) log10( (double)vector[k][i] );
                      break;
                   case 2:
                      vector[k][i] = (float) log10( (double)vector[k][i] + 1.0 );
                      break;
                   case 3:
                      vector[k][i] = (float) log( (double)vector[k][i] );
                      break;
                   case 4:
                      vector[k][i] = (float) log( (double)vector[k][i] + 1.0 );
                      break;
                   default:
                      break;
               }
           }
         }
         else if ( image->data_storage_type == VFF_TYP_COMPLEX )
         {
           points = dimension * 2;
           for ( i = 0, j = 0; i < points; i += 2, j++ ) 
           {
               vector[k][j] =  (float) cmagsq( vector[k][i], vector[k][i+1] );
               switch ( power_type )
               {
                   case 0:
                      break;
                   case 1:
                      vector[k][j] = (float) log10( (double)vector[k][j] );
                      break;
                   case 2:
                      vector[k][j] = (float) log10( (double)vector[k][j] + 1.0
);
                      break;
                   case 3:
                      vector[k][j] = (float) log( (double)vector[k][j] );
                      break;
                   case 4:
                      vector[k][j] = (float) log( (double)vector[k][j] + 1.0 );
                      break;
                   default:
                      break;
               }
           }
         }
    } /* end for */

/*
 * Call dunload to return vector to its "original" format
 */
    if ( !dunload_vector ( (char **)vector,image,(unsigned long)VFF_TYP_FLOAT, 
        numvects, dimension, process_dir ) )
    {
       fprintf (stderr,"%s:  Unable to unload vector.\n", program);
       return(0);
    }

    (void) free ( (float **) vector ); 
    return(1);

} /* end of ldpower */

/***********************************************************************
*
*  Routine Name: ldphase()
*
*          Date: 12-Jul-90 (?)
*        
*       Purpose: Phase of 1-D signal(s) 
*
*         Input: image          pointer to VIFF structure containing image 
*                               data to be processed.
*
*                phase_type     0 for radians or 1 for degrees.
*
*                procdir        process direction:  0 indicated vector 
*                               oriented processing, 1 indicates band 
*                               oriented processing.
*
*        Output: image          pointer to VIFF structure containing image 
*                               data after processing.
*
*    Written By: Ramiro Jordan 
*
* Modifications: Jeremy Worley 14-Aug-90 Incorporated into dmpp and made
*                                        static.
*                Jeremy Worley 19-Dec-90 Added case for (0,0) so that a
*                                        sqrt domain error does not occur`
*
*                Jeremy Worley 03-Nov-91 Added min_phase, which allows
*                                        us to view the radian phase in
*                                        a cycle that starts at 0 radians
*                                        or at -pi radians. 
*
***********************************************************************/

static int ldphase ( image, phase_type, min_phase,process_dir)

struct xvimage  *image;
int             phase_type,  /* phase type (radians=0 or degrees=1) */
                min_phase,   /* 0 - start at -pi, 1 start at 0 radians */
                process_dir; /* direction of processing 0=vect, 1= band */
{
    int    i, j, k, points;
    int    dunload_vector();
    int    dimension, numvects;
    float  **vector;
    char   **dload_vector();

    char *program = "ldphase";

    /* perform error checking */
    if ( !(propertype(program, image, VFF_TYP_FLOAT, FALSE)) &&
         !(propertype(program, image, VFF_TYP_COMPLEX, FALSE)) )
    {
          (void) fprintf (stderr, "\n\n%s:   ", program);
          (void) fprintf (stderr, "ldphase: data storage type must be of type, float or complex\n\n");
          return (0);
    }

/*
 * Call dload_vector for input image
 */
    if ( (vector = (float **)dload_vector(image,&numvects,&dimension,process_dir)) == NULL )
    {
       fprintf(stderr,"%s:  Unable to load vector.\n", program);
       return(0);
    }

/*
 * Perform the desired operation on the data sequence
 */
   for ( k = 0; k < numvects; k++ )
   {
        if ( image->data_storage_type == VFF_TYP_FLOAT )
        {
           for ( i = 0; i < dimension; i++ )
           {
               vector[k][i] =  0.0;
           }
         }
         else if ( image->data_storage_type == VFF_TYP_COMPLEX )
         {
           points = dimension * 2;
           for ( i = 0, j = 0; i < points; i += 2, j++ ) 
           {
               if(vector[k][i]==0.0 && vector[k][i+1]==0.0){
                  vector[k][j] = 0.0;
               }else{
                  vector[k][j] = (float) cang ( vector[k][i], vector[k][i+1] ); 
               }

               if (min_phase == 1 && vector[k][j] < 0.0){
                  vector[k][j] = vector[k][j] + XV_2PI; 
               }

               if ( phase_type == 1 ) {
                  vector[k][j] = vector[k][j] * CONSTANT;
               }
           }
         }
    } /* end for */

/*
 * Call dunload to return vector to its "original" format
 */
    if ( !dunload_vector ( (char **)vector, image,(unsigned long)VFF_TYP_FLOAT,
        numvects, dimension, process_dir ) )
    {
       fprintf (stderr,"%s:  Unable to unload vector.\n", program);
       return(0);
    }

    return(1);

} /* end of ldphase */
/* -library_code_end */
