/*****************************************************************************
 *
 * This software module was originally developed by
 *
 *   Noel Brady (Teltec DCU / ACTS-MoMuSyS)xxx (yyy / ACTS-MoMuSyS)
 *
 * and edited by
 * 
 *   Noel O'Connor (Teltec DCU / ACTS-MoMuSyS)
 *   Robert Danielsen (Telenor / ACTS-MoMuSyS)
 *   Michael Wollborn (TUH / ACTS-MoMuSys)
 *   Minhua Zhou (HHI / ACTS-MoMuSys)
 *   Michael Frater (UNSW)
 *
 * in the course of development of the MPEG-4 Video (ISO/IEC 14496-2) standard.
 * This software module is an implementation of a part of one or more MPEG-4
 * Video (ISO/IEC 14496-2) tools as specified by the MPEG-4 Video (ISO/IEC
 * 14496-2) standard.
 *
 * ISO/IEC gives users of the MPEG-4 Video (ISO/IEC 14496-2) standard free
 * license to this software module or modifications thereof for use in hardware
 * or software products claiming conformance to the MPEG-4 Video (ISO/IEC
 * 14496-2) standard.
 *
 * Those intending to use this software module in hardware or software products
 * are advised that its use may infringe existing patents. The original
 * developer of this software module and his/her company, the subsequent
 * editors and their companies, and ISO/IEC have no liability for use of this
 * software module or modifications thereof in an implementation. Copyright is
 * not released for non MPEG-4 Video (ISO/IEC 14496-2) Standard conforming
 * products.
 *
 * ACTS-MoMuSys partners retain full right to use the code for his/her own
 * purpose, assign or donate the code to a third party and to inhibit third
 * parties from using the code for non MPEG-4 Video (ISO/IEC 14496-2) Standard
 * conforming products. This copyright notice must be included in all copies or
 * derivative works.
 *
 * Copyright (c) 1996
 *
 *****************************************************************************/

/***********************************************************HeaderBegin*******
 *                                                                         
 * File:	mom_image.c
 *
 * Author:	Noel Brady (Teltec-DCU)
 * Created:	20-Feb-96
 *                                                                         
 * Description: Operations for MOMUSYS Image Structure.
 *
 * Notes: 	
 *
 * Modified:	7-2-96  Noel Brady: Access functions used.
 *		7-2-96  Noel Brady: PSNR,MSE,MAD implemented
 *		8-2-96  Noel Brady: New allocation functions
 *			i.e. Realloc..,AllocSame..
 *		20-2-96 Noel Brady/Michael Wollborn: Support Functions for
 *			generic type Image added.
 *		20-2-96 Noel Brady: Main functions are commented using header
 *		11-3-96 Noel O'Connor : Included clipping routines written
 *			by Cor Quist (KPN)
 * 		21.04.96 Robert Danielsen: Reformatted. New headers.
 *		08.05.96 Noel O'Connor: New function - ValidCoordinate()
 *    2-Dec-96 Jan De Lameillieure (HHI) : calculate SNR only in pels that belong to 
 *                                           original AND reconstructed segment mask
 *
 *		17.04.1997 Minhua (HHI) added Clean_rec_VOP
 *		07.05.97 Noel Brady added PutSubImage and modified GetSubImage
 *              03.01.98 Michael Frater: support for non-8-bit video
 *
 ***********************************************************HeaderEnd*********/

/************************    INCLUDE FILES    ********************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "momusys.h"
#include "mom_access.h"
#include "mom_image.h"
#include "vm_compos.h"

/*****
*		Image Allocation Functions 
*
*
*****/

/***********************************************************CommentBegin******
 *
 * -- AllocImage -- Allocates memory for an Image structure
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Allocates memory for an Image structure, including memory
 *	for the actual image pixels. The image pixels can be of type
 *	UChar, SInt or Float depending on the parameter	type.
 * 
 * Arguments in : 	
 *	UInt 	   size_x  : image width in pixels
 * 	UInt	   size_y  : image height in pixels
 * 	ImageType  type    : type of pixel {SHORT_TYPE, FLOAT_TYPE,
 *                                          UCHAR_TYPE}
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	A pointer to the allocated Image structure.
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	Allocates memory for the Image if possible, otherwise a 
 *	memory allocation error is signalled and the program will
 *	exit.
 *
 * See also :
 *	AllocSameImage,ReallocImage,FreeImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Image *
AllocImage(UInt size_x, UInt size_y, ImageType type)
{
  Image *image;
  
  image = (Image *) emalloc(sizeof(Image));
  
  image->version = VERSION;
  image->x = size_x;
  image->y = size_y;
  image->upperodd = 0;
  image->grid = 's';
  image->type = type;
  image->data = (ImageData *) emalloc(sizeof(ImageData)); 
  
  switch(type)
    {
    case SHORT_TYPE:
      image->data->s = (SInt *) ecalloc(size_x*size_y,sizeof(SInt));
      break;

    case FLOAT_TYPE:
      image->data->f = (Float *) ecalloc(size_x*size_y,sizeof(Float));
      break;

    case UCHAR_TYPE:
      image->data->u = (UChar *) ecalloc(size_x*size_y,sizeof(UChar));
      break;
      
    default:
      printf("Image type >>%d<< not supported\n",type);
      return(NULL);
    }
  
  image->f = image->data->s;	/* For compatibility with KHOROS */
  
  return(image);
}



/***********************************************************CommentBegin******
 *
 * -- AllocSameImage -- Allocates memory for a new Image structure
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Allocates memory for a new Image structure which is of
 * 	the same size and type as the given input image.
 * 
 * Arguments in : 	
 *	Image *image : An image which was previously allocated
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Allocates memory for the Image if possible, otherwise a 
 *	memory allocation error is signalled and the program will
 *	exit.
 *
 * See also :
 *	AllocImage, FreeImage, ReallocImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Image *
AllocSameImage(Image *image)
{
  return(AllocImage(GetImageSizeX(image),GetImageSizeY(image),
		    GetImageType(image)));
}

/***********************************************************CommentBegin******
 *
 * -- ReallocImage -- Reallocates memory for an Image structure
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Reallocates memory for an Image structure which was previously
 * 	allocated.
 * 
 * Arguments in : 	
 *	UInt size_x : the new image width
 * 	UInt size_y : the new image height
 *
 * Arguments in/out :	
 *	Image *image : the image to be reallocated memory
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	none
 *
 * Side effects :	
 *	The image structure has a new size.
 *
 * Description :	
 *	Reallocates memory if possible otherwise a memory allocation
 *	error is signalled and the program exits.
 *
 * See also :
 *	AllocImage,AllocSameImage, FreeImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
ReallocImage(Image *image, UInt size_x, UInt size_y)
{
  SInt	*ps;
  Float	*pf;
  UChar 	*pu;

  switch(GetImageType(image))
    {
    case SHORT_TYPE:
      ps = (SInt *)GetImageData(image);
      image->x = size_x;
      image->y = size_y;
      ps = (SInt *)erealloc((Char *)ps,size_x*size_y*sizeof(SInt));
      image->data->s = image->f = ps;
      break;

    case FLOAT_TYPE:
      pf = (Float *)GetImageData(image);
      image->x = size_x;
      image->y = size_y;
      pf = (Float *)erealloc((Char *)pf,size_x*size_y*sizeof(Float));
      image->data->f = pf;
      break;

    case UCHAR_TYPE:
      pu = (UChar *)GetImageData(image);
      image->x = size_x;
      image->y = size_y;
      pu = (UChar *)erealloc((Char *)pu,size_x*size_y*
			     sizeof(UChar));
      image->data->u = pu;
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image));
      return;
    }
}



/***********************************************************CommentBegin******
 *
 * -- FreeImage -- Frees up all the memory associated with an Image structure
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Frees up all the memory associated with an Image structure
 * 
 * Arguments in : 	
 *	Image *image : the Image structure to be freed
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	AllocImage,AllocSameImage,ReallocImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
FreeImage(Image *image)
{
  SInt	*ps;
  Float	*pf;
  UChar 	*pu;

  if (image == NULL) return;

  switch(GetImageType(image))
    {
    case SHORT_TYPE:
      ps = (SInt *)GetImageData(image);
      if(ps != NULL) free((Char *)ps);
      free((Char *) image->data);
      free((Char *)image);
      break;

    case FLOAT_TYPE:
      pf = (Float *)GetImageData(image);
      if(pf != NULL) free((Char *)pf);
      free((Char *) image->data);
      free((Char *)image);
      break;

    case UCHAR_TYPE:
      pu = (UChar *)GetImageData(image);
      if(pu != NULL) free((Char *)pu);
      free((Char *) image->data);
      free((Char *)image);
      break;
      
    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- Alloc{xxx} -- Allocation functions to be discontinued
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Allocation functions to be discontinued
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

ImageI *
AllocImageI(UInt size_x, UInt size_y)
{
  return(AllocImage(size_x, size_y, SHORT_TYPE));
}

ImageF *
AllocImageF(UInt size_x, UInt size_y)
{
  return(AllocImage(size_x, size_y, FLOAT_TYPE));
}

ImageI *
AllocSameImageI(ImageI *image)
{
  return AllocImage(GetImageSizeX(image),GetImageSizeY(image),SHORT_TYPE);
}

ImageF *
AllocSameImageF(ImageF *image)
{
  return AllocImage(GetImageSizeX(image),GetImageSizeY(image),FLOAT_TYPE);
}


/***********************************************************CommentBegin******
 *
 * -- ReallocImage{x} -- Realloc functions for ImageI and ImageF
 *
 *	Void ReallocImageI(ImageI *image, UInt size_x, UInt size_y)
 *	Void ReallocImageF(ImageF *image, UInt size_x, UInt size_y)
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Realloc functions for ImageI and ImageF
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
ReallocImageI(ImageI *image, UInt size_x, UInt size_y)
{
  ReallocImage(image,size_x,size_y);
}

Void
ReallocImageF(ImageF *image, UInt size_x, UInt size_y)
{
  ReallocImage(image,size_x,size_y);
}



/***********************************************************CommentBegin******
 *
 * -- FreeImage{x} -- Freeing functions for ImageI and ImageF
 *
 *	Void FreeImageI(ImageI *image)
 *	Void FreeImageF(ImageF *image)
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Freeing functions for ImageI and ImageF
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
FreeImageI(ImageI *image)
{
  FreeImage(image);
}

Void
FreeImageF(ImageF *image)
{
  FreeImage(image);
}



/***********************************************************CommentBegin******
 *
 * -- CopyImage -- Copies the pixel values of one image to another image
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Copies the pixel values of one image to another image.
 * 
 * Arguments in : 	
 *	Image *image_in : the source image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *image_out : the destination image
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, the function exits with an error if
 *	the source and destination images are not of the same allocated
 *	size and type.
 *
 *	Also, the images of type UChar are not yet supported.
 *
 * See also :
 *	CastImageItoF, CastImageFtoI
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
CopyImage(Image *image_in, Image *image_out)
{
  switch(GetImageType(image_out))
    {
    case SHORT_TYPE:
      CopyImageI(image_in,image_out);
      break;

    case FLOAT_TYPE:
      CopyImageF(image_in,image_out);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image_out));
      return;
    }
}



/***********************************************************CommentBegin******
 *
 * -- CopyImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :	
 *	    26.05.97 Minhua Zhou: allowed to copy images with different
 *                                size; 
 *	
 *
 ***********************************************************CommentEnd********/

Void
CopyImageI(ImageI *image_in, ImageI *image_out)
{
  SInt 	*p_in  = GetImageIData(image_in),
    *p_out = GetImageIData(image_out),
    *p_end;

  UInt	sx_in  = GetImageISizeX(image_in),
    sx_out = GetImageISizeX(image_out),
    sy_in  = GetImageISizeY(image_in),
    sy_out = GetImageISizeY(image_out),
    sxy_in = GetImageISize(image_in);

  if (GetImageType(image_in) != SHORT_TYPE) {
    runerr("ERROR(CopyImageI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(CopyImageI): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      Char message[101];
      sprintf(message,
             "ERROR(CopyImageI): images are different sizes %ldx%ld - %ldx%ld",
             sx_in,sy_in,sx_out,sy_out);
      runerr(message,1);
        }

  p_end = p_in + sxy_in;

  while (p_in != p_end) 
    {
      *p_out = *p_in;
      p_in++;
      p_out++;
    }

}



/***********************************************************CommentBegin******
 *
 * -- CopyImageF -- Copies image_in to image_out
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Copies image_in to image_out
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	An error is signalled if the images do not have the same dimensions.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
CopyImageF(ImageF *image_in, ImageF *image_out)
{
  Float *p_in  = GetImageFData(image_in),
    *p_out = GetImageFData(image_out),
    *p_end;

  UInt    sx_in = GetImageFSizeX(image_in),
    sx_out = GetImageFSizeX(image_out),
    sy_in = GetImageFSizeY(image_in),
    sy_out = GetImageFSizeY(image_out),
    sxy_in = GetImageFSize(image_in);


  if (GetImageType(image_in) != FLOAT_TYPE) {
    runerr("ERROR(CopyImageF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(CopyImageF): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      runerr("ERROR(CopyImageF): images are different sizes",1);
    }

  p_end = p_in + sxy_in;

  while (p_in != p_end) 
    {
      *p_out = *p_in;
      p_in++;
      p_out++;
    }

}


/***********************************************************CommentBegin******
 *
 * -- CastImageItoF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
CastImageItoF(ImageI *image_in, ImageF *image_out)
{
  SInt 	*p_in  = GetImageIData(image_in),
    *p_end;

  Float	*p_out = GetImageFData(image_out);

  UInt	sx_in  = GetImageISizeX(image_in),
    sx_out = GetImageFSizeX(image_out),
    sy_in  = GetImageISizeY(image_in),
    sy_out = GetImageFSizeY(image_out),
    sxy_in = GetImageISize(image_in);


  if (GetImageType(image_in) != SHORT_TYPE) {
    runerr("ERROR(CastImageItoF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(CastImageItoF): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      runerr("ERROR(CastImageItoF): images are different sizes",1);
    }

  p_end = p_in + sxy_in;

  while (p_in != p_end) 
    {
      *p_out = (Float)(*p_in);
      p_in++;
      p_out++;
    }

}


/***********************************************************CommentBegin******
 *
 * -- CastImageFtoI -- Casts ImageF to ImageI
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Casts image_in of type ImageF to image_out of type ImageI
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	An error is signalled if the images do not have the same dimensions.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
CastImageFtoI(ImageF *image_in, ImageI *image_out)
{
  Float	*p_in  = GetImageFData(image_in),
    *p_end;

  SInt 	*p_out = GetImageIData(image_out);

  UInt	sx_in  = GetImageFSizeX(image_in),
    sx_out = GetImageISizeX(image_out),
    sy_in  = GetImageFSizeY(image_in),
    sy_out = GetImageISizeY(image_out),
    sxy_in = GetImageFSize(image_in);

  if (GetImageType(image_in) != FLOAT_TYPE) {
    runerr("ERROR(CastImageFtoI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(CastImageFtoI): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      runerr("ERROR(CastImageFtoI): images are different sizes",1);
    }

  p_end = p_in + sxy_in;

  while (p_in != p_end) 
    {
      *p_out = (SInt)(*p_in);
      p_in++;
      p_out++;
    }

}


/*****	
*	
*			Arithmetic Image functions 
*
*****/


/***********************************************************CommentBegin******
 *
 * -- SetConstantImage -- Sets each pixel in the Image to val
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Sets each pixel in the Image to val.
 * 
 * Arguments in : 	
 *	Float	val : the pixel value to be used
 *
 * Arguments in/out :	
 *	Image *image: an image structure
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	Note val is a Float. If image is not a Float image then
 *	val is cast to the appropriate type i.e. SInt, unsigned
 *	Char.
 *
 * Description :	
 *	Images of type UChar are not yet supported.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SetConstantImage(Image *image, Float val)
{
  switch(GetImageType(image))
    {
    case SHORT_TYPE:
      SetConstantImageI(image,(SInt)val);
      break;

    case FLOAT_TYPE:
      SetConstantImageF(image,val);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image));
      return;
    }
}



/***********************************************************CommentBegin******
 *
 * -- SetConstantImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SetConstantImageI(ImageI *image, SInt val)
{
  SInt	*p  = GetImageIData(image),
    *p_end;

  UInt	sxy = GetImageISize(image);

  if (GetImageType(image) != SHORT_TYPE) {
    runerr("ERROR(SetConstantImageI): image not of expected type",1);
  }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = val;
      p++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- SetConstantImageF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SetConstantImageF(ImageF *image, Float val)
{
  Float		*p  = GetImageFData(image),
    *p_end;

  UInt	sxy = GetImageFSize(image);

  if (GetImageType(image) != FLOAT_TYPE) {
    runerr("ERROR(SetConstantImageF): image not of expected type",1);
  }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = val;
      p++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AbsImage -- Copies pixel absolute values of one image to another image
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Copies the pixel absolute values of one image to another image.
 * 
 * Arguments in : 	
 *	Image *image_in : the source image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *image_out : the destination image
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, the function exits with an error if
 *	the source and destination images are not of the same allocated
 *	size and type.
 *
 *	Also, the images of type UChar are not yet supported.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/
 
Void
AbsImage(Image *image_in, Image *image_out)
{
  switch(GetImageType(image_in))
    {
    case SHORT_TYPE:
      AbsImageI(image_in,image_out);
      break;

    case FLOAT_TYPE:
      AbsImageF(image_in,image_out);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image_in));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- AbsImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AbsImageI(ImageI *image_in, ImageI *image_out)
{
  SInt	*p_in  = GetImageIData(image_in),
    *p_out = GetImageIData(image_out),
    *p_end;

  UInt	sx_in  = GetImageISizeX(image_in),
    sx_out = GetImageISizeX(image_out),
    sy_in  = GetImageISizeY(image_in),
    sy_out = GetImageISizeY(image_out),
    sxy_in = GetImageISize(image_in);

  if (GetImageType(image_in) != SHORT_TYPE) {
    runerr("ERROR(AbsImageI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(AbsImageI): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      runerr("ERROR(AbsImageI): images are different sizes",1);
    }

  p_end = p_in + sxy_in;

  while (p_in != p_end)
    {
      *p_out = abs(*p_in);
      p_in++;
      p_out++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AbsImageF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AbsImageF(ImageF *image_in, ImageF *image_out)
{
  Float		*p_in  = GetImageFData(image_in),
    *p_out = GetImageFData(image_out),
    *p_end;

  UInt	sx_in  = GetImageFSizeX(image_in),
    sx_out = GetImageFSizeX(image_out),
    sy_in  = GetImageFSizeY(image_in),
    sy_out = GetImageFSizeY(image_out),
    sxy_in = GetImageFSize(image_in);

  if (GetImageType(image_in) != FLOAT_TYPE) {
    runerr("ERROR(AbsImageF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(AbsImageF): image not of expected type",1);
  }

  if ((sx_in != sx_out) || (sy_in != sy_out)) 
    {
      runerr("ERROR(AbsImageF): images are different sizes",1);
    }

  p_end = p_in + sxy_in;

  while (p_in != p_end)
    {
      *p_out = fabs(*p_in);
      p_in++;
      p_out++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AddImage -- Adds two images
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Adds two images and stores the result in the third
 * 	i.e. image_out = image_in1 + image_in2
 * 
 * Arguments in : 	
 *	Image *image_in1 : the first source image
 *	Image *image_in12 : the second source image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *image_out : the image storing the addition result
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, the function exits with an error if
 *	the source and destination images are not of the same allocated
 *	size and type.
 *
 *	Also, the images of type UChar are not yet supported.
 *
 * See also :
 *	SubImage, AbsSubImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AddImage(Image *image_in1, Image *image_in2, Image *image_out)
{
  switch(GetImageType(image_in1))
    {
    case SHORT_TYPE:
      AddImageI(image_in1,image_in2,image_out);
      break;

    case FLOAT_TYPE:
      AddImageF(image_in1,image_in2,image_out);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image_in1));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- AddImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AddImageI(ImageI *image_in1, ImageI *image_in2, ImageI *image_out)
{
  SInt	*p  = GetImageIData(image_out),
    *p1 = GetImageIData(image_in1),
    *p2 = GetImageIData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageISizeX(image_in1),
    sx_in2 = GetImageISizeX(image_in2),
    sx_out = GetImageISizeX(image_out),
    sy_in1 = GetImageISizeY(image_in1),
    sy_in2 = GetImageISizeY(image_in2),
    sy_out = GetImageISizeY(image_out),
    sxy    = GetImageISize(image_out);

  if (GetImageType(image_in1) != SHORT_TYPE) {
    runerr("ERROR(AddImageI): image not of expected type",1);
  }

  if (GetImageType(image_in2) != SHORT_TYPE) {
    runerr("ERROR(AddImageI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(AddImageI): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(AddImageI): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(AddImageI): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = *p1 + *p2;
      p++;
      p1++;
      p2++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AddImageF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AddImageF(ImageF *image_in1, ImageF *image_in2, ImageF *image_out)
{
  Float	*p  = GetImageFData(image_out),
    *p1 = GetImageFData(image_in1),
    *p2 = GetImageFData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageFSizeX(image_in1),
    sx_in2 = GetImageFSizeX(image_in2),
    sx_out = GetImageFSizeX(image_out),
    sy_in1 = GetImageFSizeY(image_in1),
    sy_in2 = GetImageFSizeY(image_in2),
    sy_out = GetImageFSizeY(image_out),
    sxy    = GetImageFSize(image_out);

  if (GetImageType(image_in1) != FLOAT_TYPE) {
    runerr("ERROR(AddImageF): image not of expected type",1);
  }

  if (GetImageType(image_in2) != FLOAT_TYPE) {
    runerr("ERROR(AddImageF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(AddImageF): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(AddImageF): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(AddImageF): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = *p1 + *p2;
      p++;
      p1++;
      p2++;
    }

}
 
/***********************************************************CommentBegin******
 *
 * -- SubImage -- Subtracts two images
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Subtracts two images and stores the result in the third
 * 	i.e. image_out = image_in1 - image_in2
 * 
 * Arguments in : 	
 *	Image *image_in1 : the first source image
 *	Image *image_in12 : the second source image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *image_out : the image storing the subtraction result
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, the function exits with an error if
 *	the source and destination images are not of the same allocated
 *	size and type.
 *
 *	Also, the images of type UChar are not yet supported.
 *
 * See also :
 *	AddImage, AbsSubImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SubImage(Image *image_in1, Image *image_in2, Image *image_out)
{
  switch(GetImageType(image_in1))
    {
    case SHORT_TYPE:
      SubImageI(image_in1,image_in2,image_out);
      break;

    case FLOAT_TYPE:
      SubImageF(image_in1,image_in2,image_out);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image_in1));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- SubImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SubImageI(ImageI *image_in1, ImageI *image_in2, ImageI *image_out)
{
  SInt		*p  = GetImageIData(image_out),
    *p1 = GetImageIData(image_in1),
    *p2 = GetImageIData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageISizeX(image_in1),
    sx_in2 = GetImageISizeX(image_in2),
    sx_out = GetImageISizeX(image_out),
    sy_in1 = GetImageISizeY(image_in1),
    sy_in2 = GetImageISizeY(image_in2),
    sy_out = GetImageISizeY(image_out),
    sxy    = GetImageISize(image_out);

  if (GetImageType(image_in1) != SHORT_TYPE) {
    runerr("ERROR(SubImageI): image not of expected type",1);
  }

  if (GetImageType(image_in2) != SHORT_TYPE) {
    runerr("ERROR(SubImageI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(SubImageI): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(SubImageI): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(SubImageI): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = *p1 - *p2;
      p++;
      p1++;
      p2++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- SubImageF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
SubImageF(ImageF *image_in1, ImageF *image_in2, ImageF *image_out)
{
  Float	*p  = GetImageFData(image_out),
    *p1 = GetImageFData(image_in1),
    *p2 = GetImageFData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageFSizeX(image_in1),
    sx_in2 = GetImageFSizeX(image_in2),
    sx_out = GetImageFSizeX(image_out),							
    sy_in1 = GetImageFSizeY(image_in1),
    sy_in2 = GetImageFSizeY(image_in2),
    sy_out = GetImageFSizeY(image_out),
    sxy    = GetImageFSize(image_out);

  if (GetImageType(image_in1) != FLOAT_TYPE) {
    runerr("ERROR(SubImageF): image not of expected type",1);
  }

  if (GetImageType(image_in2) != FLOAT_TYPE) {
    runerr("ERROR(SubImageF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(SubImageF): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(SubImageF): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(SubImageF): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = *p1 - *p2;
      p++;
      p1++;
      p2++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AbsSubImage -- Subtracts two images
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Subtracts two images and stores the absolute result in the third
 * 	i.e. image_out = abs(image_in1 - image_in2)
 * 
 * Arguments in : 	
 *	Image *image_in1 : the first source image
 *	Image *image_in12 : the second source image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *image_out : the image storing the subtraction result
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, the function exits with an error if
 *	the source and destination images are not of the same allocated
 *	size and type.
 *
 *	Also, the images of type UChar are not yet supported.
 *
 * See also :
 *	AddImage, SubImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AbsSubImage(Image *image_in1, Image *image_in2, Image *image_out)
{
  switch(GetImageType(image_in1))
    {
    case SHORT_TYPE:
      AbsSubImageI(image_in1,image_in2,image_out);
      break;

    case FLOAT_TYPE:
      AbsSubImageF(image_in1,image_in2,image_out);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image_in1));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- AbsSubImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AbsSubImageI(ImageI *image_in1, ImageI *image_in2, ImageI *image_out)
{
  SInt	*p  = GetImageIData(image_out),
    *p1 = GetImageIData(image_in1),
    *p2 = GetImageIData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageISizeX(image_in1),
    sx_in2 = GetImageISizeX(image_in2),
    sx_out = GetImageISizeX(image_out),
    sy_in1 = GetImageISizeY(image_in1),
    sy_in2 = GetImageISizeY(image_in2),
    sy_out = GetImageISizeY(image_out),
    sxy    = GetImageISize(image_out);

  if (GetImageType(image_in1) != SHORT_TYPE) {
    runerr("ERROR(AbsSubImageI): image not of expected type",1);
  }

  if (GetImageType(image_in2) != SHORT_TYPE) {
    runerr("ERROR(AbsSubImageI): image not of expected type",1);
  }

  if (GetImageType(image_out) != SHORT_TYPE) {
    runerr("ERROR(AbsSubImageI): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(AbsSubImageI): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(AbsSubImageI): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = abs(*p1 - *p2);
      p++;
      p1++;
      p2++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- AbsSubImageF -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
AbsSubImageF(ImageF *image_in1, ImageF *image_in2, ImageF *image_out)
{
  Float	*p  = GetImageFData(image_out),
    *p1 = GetImageFData(image_in1),
    *p2 = GetImageFData(image_in2),
    *p_end;

  UInt	sx_in1 = GetImageFSizeX(image_in1),
    sx_in2 = GetImageFSizeX(image_in2),
    sx_out = GetImageFSizeX(image_out),
    sy_in1 = GetImageFSizeY(image_in1),
    sy_in2 = GetImageFSizeY(image_in2),
    sy_out = GetImageFSizeY(image_out),
    sxy    = GetImageFSize(image_out);

  if (GetImageType(image_in1) != FLOAT_TYPE) {
    runerr("ERROR(AbsSubImageF): image not of expected type",1);
  }

  if (GetImageType(image_in2) != FLOAT_TYPE) {
    runerr("ERROR(AbsSubImageF): image not of expected type",1);
  }

  if (GetImageType(image_out) != FLOAT_TYPE) {
    runerr("ERROR(AbsSubImageF): image not of expected type",1);
  }

  if ((sx_in1 != sx_in2) || (sy_in1 != sy_in2)) 
    {
      runerr("ERROR(AbsSubImageF): images are different sizes",1);
    }
  if ((sx_in1 != sx_out) || (sy_in1 != sy_out)) 
    {
      runerr("ERROR(AbsSubImageF): images are different sizes",1);
    }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = fabs(*p1 - *p2);
      p++;
      p1++;
      p2++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- ClipImage -- To Clip an image within 0 and 255
 *
 * Author :		
 *	Cor Quist, KPN Research
 *
 * Created :		
 *	March 11, 1996
 *
 * Purpose :		
 *	To Clip an image within 0 and 255
 * 
 * Arguments in : 	
 *	Image *image
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	Image *image
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Clips an image between 0 and 255
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/
 
Void
ClipImage(Image *image, Int maxval)
{
  switch(GetImageType(image))
    {
    case SHORT_TYPE:
      ClipImageI(image, maxval);
      break;

    case FLOAT_TYPE:
      ClipImageF(image, maxval);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image));
      return;
    }
}

/***********************************************************CommentBegin******
 *
 * -- ClipImageI -- 
 *
 * Author :		
 *	Cor Quist, KPN Research
 *
 * Created :		
 *	March 11, 1996
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
ClipImageI(ImageI *image, Int maxval)
{
  SInt	*p  = GetImageIData(image),
        *p_end;

  ULInt	sxy = GetImageISize(image);

  if (GetImageType(image) != SHORT_TYPE) {
    runerr("ERROR(ClipImageI): image not of expected type",1);
  }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = CLIP(*p, 0, maxval);
      p++;
    }

}

/***********************************************************CommentBegin******
 *
 * -- ClipImageF -- 
 *
 * Author :		
 *	Cor Quist, KPN Research
 *
 * Created :		
 *	March 11, 1996
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
ClipImageF(ImageF *image, Int maxval)
{
  Float	       *p  = GetImageFData(image),
               *p_end;

  ULInt	sxy = GetImageFSize(image);

  if (GetImageType(image) != FLOAT_TYPE) {
    runerr("ERROR(ClipImageF): image not of expected type",1);
  }

  p_end = p + sxy;

  while (p != p_end)
    {
      *p = CLIP(*p, 0, maxval);
      p++;
    }

}

/*  ***** Image Quality Functions **** */

/***********************************************************CommentBegin******
 *
 * -- PSNRImage -- Computes the PSNR of an image with respect to another
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Computes the PSNR of an image with respect to another
 * 
 * Arguments in : 	
 *	Image *image1 : the first image
 *	Image *image2 : the second image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	the peak signal-to-noise ratio
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, only images of type SInt are supported.
 *
 * See also :
 *	MSEImage, MADImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
PSNRImage(Image *image1, Image *image2)
{
  Float	mse, psnr;

  mse = MSEImage(image1,image2);
  psnr = 20*log10(((Float)255.0)/sqrt(mse));

  return psnr;
}



/***********************************************************CommentBegin******
 *
 * -- MaskedPSNRImage -- Computes the PSNR of an masked image
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Computes the PSNR of a masked image
 * 
 * Arguments in : 	
 *	Image *image1 : the first image
 *	Image *image2 : the second image
 *	Image *mask : the mask
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	the peak signal-to-noise ratio
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, only images of type SInt are supported.
 *
 * See also :
 *
 *
 * Modified :		
 *	2-Dec-96 Jan De Lameillieure (HHI) : calculate SNR only in pels that belong to 
 *                                           original AND reconstructed segment mask
 *
 *      23-Jul-97 Ulrike Pestel-Schiller(UH):added warning for empty mask
 ***********************************************************CommentEnd********/

Float
MaskedPSNRImage(Image *image1, Image *image2, Image *mask1, Image *mask2)
{
  Float	mse, psnr;
	Image 	*im1 = AllocSameImage(image1),
					*im2 = AllocSameImage(image2);
	Int 	n;

	MaskImage(image1, mask1, mask2, im1);
	n = MaskImage(image2, mask1, mask2, im2);
	if(n==0)
	  printf("*********WARNING MaskedPSNRImage: empty mask*************\n"); 

  mse = MSEImage(im1,im2) * GetImageSize(image1) / n;
  psnr = 20*log10(((Float)255.0)/sqrt(mse));

 	FreeImage(im1);
	FreeImage(im2);

  return psnr;
}




/***********************************************************CommentBegin******
 *
 * -- MSEImage -- Computes the MSE of an image with respect to another
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Computes the MSE of an image with respect to another
 * 
 * Arguments in : 	
 *	Image *image1 : the first image
 *	Image *image2 : the second image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	the mean square error
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, only images of type SInt are supported.
 *
 * See also :
 *	PSNRImage, MADImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
MSEImage(Image *image1, Image *image2)
{
  Float mse;
 	
  switch(GetImageType(image1))
    {
    case SHORT_TYPE:
      mse = MSEImageI(image1,image2);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image1));
      return 0.0;
    }
  return mse;
}

/***********************************************************CommentBegin******
 *
 * -- MSEImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
MSEImageI(ImageI *image1, ImageI *image2)
{
  ImageI 	*image3;

  SInt	*p, *p_end;

  UInt 	s, sum=0;

  Float	mse;

  if (GetImageType(image1) != SHORT_TYPE) {
    runerr("ERROR(MSEImageI): image not of expected type",1);
  }

  if (GetImageType(image2) != SHORT_TYPE) {
    runerr("ERROR(MSEImageI): image not of expected type",1);
  }

  image3 = AllocSameImageI(image1);

  SubImageI(image1,image2,image3);

  p = GetImageIData(image3);
  s = GetImageISize(image3);

  p_end = p + s;

  while (p != p_end) 
    {
      sum = sum + (*p) * (*p);
      p++;
    }
	
  mse = ((Float) sum)/s;

  FreeImageI(image3);

  return mse;
}

/***********************************************************CommentBegin******
 *
 * -- MADImage -- Computes the MAD of an image with respect to another
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Computes the MAD of an image with respect to another
 * 
 * Arguments in : 	
 *	Image *image1 : the first image
 *	Image *image2 : the second image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	the mean absolute difference
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Currently, only images of type SInt are supported.
 *
 * See also :
 *	MSEImage, PSNRImage
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
MADImage(Image *image1, Image *image2)
{
  Float mad;
    
  switch(GetImageType(image1))
    {
    case SHORT_TYPE:
      mad = MADImageI(image1,image2);
      break;

    default:
      printf("Image type >>%d<< not supported\n",GetImageType(image1));
      return 0.0;
    }
  return mad;
}

/***********************************************************CommentBegin******
 *
 * -- MADImageI -- 
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
MADImageI(ImageI *image1, ImageI *image2)
{
  ImageI 	*image3;

  SInt	*p, *p_end;

  UInt 	s, sum=0;

  Float	mad;

  if (GetImageType(image1) != SHORT_TYPE) {
    runerr("ERROR(MADImageI): image not of expected type",1);
  }

  if (GetImageType(image2) != SHORT_TYPE) {
    runerr("ERROR(MADImageI): image not of expected type",1);
  }

  image3 = AllocSameImageI(image1);

  AbsSubImageI(image1,image2,image3);

  p = GetImageIData(image3);
  s = GetImageISize(image3);

  p_end = p + s;

  while (p != p_end) 
    {
      sum = sum + (*p);
      p++;
    }
	
  mad = ((Float) sum)/s;

  FreeImageI(image3);

  return mad;
}

/***********************************************************CommentBegin******
 *
 * -- emalloc -- Memory allocation with error handling
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Memory allocation with error handling.
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	WARNING: There is a problem when you want to allocate more memory
 *	than size can handle.  Malloc gets as argument an unsigned.  It poses
 *	a problem when sizeof(unsigned) < sizeof(Char *)...
 *
 *	Code adapted from XLIM
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Char *
emalloc(Int n)
{
  Char 	*p;
	
  if(n <= 0) runerr("ERROR(emalloc): Wrong size: <= 0.",1);

  p = (Char *) malloc((UInt)n);
  if(p == NULL)
    runerr("ERROR(emalloc): memory allocation error.",1);

  return p;
}

/***********************************************************CommentBegin******
 *
 * -- ecalloc -- Memory allocation with error handling
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Memory allocation with error handling.
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	WARNING: There is a problem when you want to allocate more memory
 *	than size can handle.  Malloc gets as argument an unsigned.  It poses
 *	a problem when sizeof(unsigned) < sizeof(Char *)...
 *
 *	Code adapted from XLIM
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Char *
ecalloc(Int n, Int s)
{
  Char 	*p;
	
  if(s <= 0 || n <= 0) runerr("ERROR(ecalloc): Wrong size: <= 0.",1);
  p = (Char *) calloc((UInt)n,(UInt)s);
  if(p == NULL)
    runerr("ERROR(ecalloc): memory allocation error.",1);
	
  return p;
}

/***********************************************************CommentBegin******
 *
 * -- erealloc -- Memory re-allocation with error handling.
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Memory re-allocation with error handling
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	WARNING: There is a problem when you want to allocate more memory
 *	than size can handle.  Malloc gets as argument an unsigned.  It poses
 *	a problem when sizeof(unsigned) < sizeof(Char *)...
 *
 *	Code adapted from XLIM
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Char *
erealloc(Char *p, Int n)
{	
  if (n <= 0) runerr("ERROR(erealloc): Wrong size: <= 0.",1);
  p = (Char *) realloc(p,(UInt)n);
  if(p == NULL)
    runerr("ERROR(erealloc): memory allocation error.",1);
	
  return p;
}

/***********************************************************CommentBegin******
 *
 * -- runerr -- Run-time error handling function
 *
 * Author :		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Run-time error handling function
 * 
 * Arguments in : 	
 *	
 *
 * Arguments in/out :	
 *	
 *
 * Arguments out :	
 *	
 *
 * Return values :	
 *	
 *
 * Side effects :	
 *	
 *
 * Description :	
 *	Code adapted from XLIM
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
runerr(Char *error_msg, Int exit_code)
{
  printf("\n\n%s\n",error_msg);
  exit(exit_code);
}

/*******************        INCLUDE FILES       ******************************/

#include <stdio.h>
#include <math.h>
#include "momusys.h"
#include "mom_access.h"
/* #include "alp_code.h" */	/* commented by JDL 11-Sep-96 */


/***********************************************************CommentBegin******
 *
 * -- MultConstantImage -- The input image is multiplied with a constant "val"
 *
 * Author : Michael Wollborn (TUH)		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	The input image is multiplied with a constant "val".
 * 
 * Arguments in : 	
 *	Image *im_in	: input image
 *	Float val	: value to multiply image with [in]
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *im_out	: output image
 *
 * Return values :	
 *	none
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
MultConstantImage(Image *im_in, Image *im_out, Float val)
{
  Char *function_name = "MultConstantImage";
  
  Int i;

  /*****
   *
   *	Check image dimensions and type
   *
   *****/
  if(COMPARE_IMAGE_SIZE_EQL(im_in,im_out))
    {
      ERROR_IMAGE_SIZE_EQL(im_in,im_out);
      return;
    }

  if(COMPARE_IMAGE_TYPE_EQL(im_in,im_out))
    {
      ERROR_IMAGE_TYPE_EQL(im_in,im_out);
      return;
    }

  /*****
   *
   *	Work on correct image type
   *
   *****/
  switch(GetImageType(im_in))
    {
    case UCHAR_TYPE:
      {
	UChar *pim_in,*pim_out;
      
	/*****
	 *
	 *	Check if range of "UChar" is exceeded
	 *
	 *****/
	if(GetMaxValueImage(im_in)*val >= MAX_UCHAR_VAL)
	  {
	    ERROR("Result too large for unsigned char (UChar)");
	    return;
	  }
  
	/*****
	 *
	 *	Scan image and multiply each pel with "val"
	 *
	 *****/
	pim_in  = (UChar *)GetImageData(im_in);
	pim_out = (UChar *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  *pim_out = (UChar) ROUND((*pim_in) * val);
	break;
      }

    case SHORT_TYPE:
      {
	SInt *pim_in,*pim_out;
      
	/*****
	 *
	 *	Check if range of "SInt" is exceeded
	 *
	 *****/
	if(GetMaxValueImage(im_in)*val >= MAX_SHORT_VAL)
	  {
	    ERROR("Result too large for short Int (SInt)");
	    return;
	  }
  
	/*****
	 *
	 *	Scan image and multiply each pel with "val"
	 *
	 *****/
	pim_in  = (SInt *)GetImageData(im_in);
	pim_out = (SInt *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  *pim_out = (SInt) ROUND((*pim_in) * val);
	break;
      }

    case FLOAT_TYPE:
      {
	Float *pim_in,*pim_out;
      
	/*****
	 *
	 *	Check if range of "Float" is exceeded
	 *
	 *****/
	if(GetMaxValueImage(im_in)*val >= MAX_FLOAT_VAL)
	  {
	    ERROR("Result too large for float (Float)");
	    return;
	  }
  
	/*****
	 *
	 *	Scan image and multiply each pel with "val"
	 *
	 *****/
	pim_in  = (Float *)GetImageData(im_in);
	pim_out = (Float *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  *pim_out = (Float) (*pim_in) * val;
	break;
      }
      
    default:
      ERROR("Image Type not supported");
      return;
    }
  
  return;
}

/***********************************************************CommentBegin******
 *
 * --  GetMaxValueImage -- Looks for the maximum value in the image
 *
 * Author : Michael Wollborn (TUH)		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Looks for the maximum value in the image
 * 
 * Arguments in : 	
 *	Image *im_in	: input image
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	none
 *
 * Return values :	
 *	Maximum value in the image
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Float
GetMaxValueImage(Image *im_in)
{
  Char *function_name = "GetMaxValueImage";
  
  Int i;

  /*****
   *
   *	Work on correct image type
   *
   *****/
  switch(GetImageType(im_in))
    {
    case UCHAR_TYPE:
      {
	UChar *pim_in,maxval;

	pim_in  = (UChar *)GetImageData(im_in);
      
	for(i=GetImageSize(im_in),maxval=0; i; i--,pim_in++)
	  if(*pim_in > maxval)
	    maxval = *pim_in;
  
	return((Float)maxval);
      }

    case SHORT_TYPE:
      {
	SInt *pim_in,maxval;

	pim_in  = (SInt *)GetImageData(im_in);
      
	for(i=GetImageSize(im_in),maxval=0; i; i--,pim_in++)
	  if(*pim_in > maxval)
	    maxval = *pim_in;
  
	return((Float)maxval);
      }

    case FLOAT_TYPE:
      {
	Float *pim_in,maxval;

	pim_in  = (Float *)GetImageData(im_in);
      
	for(i=GetImageSize(im_in),maxval=0; i; i--,pim_in++)
	  if(*pim_in > maxval)
	    maxval = *pim_in;
  
	return((Float)maxval);
      }

    default:
      ERROR("Image Type not supported");
      return(0);
    }
}

/***********************************************************CommentBegin******
 *
 * -- GetSubImage -- Extracts a subimage from an image
 *
 * Author : Michael Wollborn (TUH)		
 *     
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Extracts a subimage from an image.
 * 
 * Arguments in : 	
 *	Image *im_in	: input image
 *	Int xpos	: left column of subimage
 *	Int ypos	: top row of subimage
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *im_out	: output image
 *
 * Return values :	
 *	none
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	A subimage is extracted from the input image. The position
 *	of the top-left pel of the subimage to extract is given
 *      by the coordinates xpos,ypos.
 *
 * See also :
 *	
 *
 * Modified :		 
 *		07.05.97 Noel Brady: Modified in order to deal with cases when the
 *													the sub-image requires pixels outside the 
 *													imput image.
 *	
 *
 ***********************************************************CommentEnd********/

Void
GetSubImage(Image *im_in, Image *im_out, Int xpos, Int ypos)
{
  Char *function_name = "GetSubImage";
  
  Int x,y,xdim_in,ydim_in,xdim_out,ydim_out,xdim_off,xdim_off1,xdim_off2;

	Int xs,ys,xe,ye;

  /*****
   *
   *	Check image dimensions and type (if im_in is >= im_out)
   *
   *****/
 

  if(COMPARE_IMAGE_TYPE_EQL(im_in,im_out))
    {
      ERROR_IMAGE_TYPE_EQL(im_in,im_out);
      return;
    }

  /*****
   *
   *	Get image dimensions
   *
   *****/
  xdim_in  = GetImageSizeX(im_in);
  ydim_in  = GetImageSizeY(im_in);

  xdim_out = GetImageSizeX(im_out);
  ydim_out = GetImageSizeY(im_out);
  
  xdim_off = xdim_in - xdim_out;

  /*****
   *
   *	Work on correct image type
   *
   *****/
  switch(GetImageType(im_in))
    {
    case UCHAR_TYPE:
      {
	UChar *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *
	 *****/
	pim_in  = (UChar *)GetImageData(im_in) + MAX(0,ypos)*xdim_in 
											+ MAX(0,xpos);
	xs = MAX(0,xpos) - xpos;
	ys = MAX(0,ypos) - ypos;
	xe = MIN(xpos+xdim_out,xdim_in)-xpos;
	ye = MIN(ypos+ydim_out,ydim_in)-ypos;
	xdim_off1 = xdim_in - (xe-xs);
	xdim_off2 = xdim_out - (xe-xs);

	pim_out = (UChar *)GetImageData(im_out) + ys*xdim_out + xs;

	for(y=ys; y<ye; y++,pim_in+=xdim_off1,pim_out+=xdim_off2)
	  for(x=xs; x<xe; x++,pim_in++,pim_out++)
	    *pim_out = *pim_in;
	break;
      }

    case SHORT_TYPE:
      {
	SInt *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *

 		*****/
	pim_in  = (SInt *)GetImageData(im_in) + MAX(0,ypos)*xdim_in 
											+ MAX(0,xpos);
	xs = MAX(0,xpos) - xpos;
	ys = MAX(0,ypos) - ypos;
	xe = MIN(xpos+xdim_out,xdim_in)-xpos;
	ye = MIN(ypos+ydim_out,ydim_in)-ypos;
	xdim_off1 = xdim_in - (xe-xs);
	xdim_off2 = xdim_out - (xe-xs);

	pim_out = (SInt *)GetImageData(im_out) + ys*xdim_out + xs;

	for(y=ys; y<ye; y++,pim_in+=xdim_off1,pim_out+=xdim_off2)
	  for(x=xs; x<xe; x++,pim_in++,pim_out++)
	    *pim_out = *pim_in;
	break;
      }

    case FLOAT_TYPE:
      {
	Float *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *
	 *****/
	pim_in  = (Float *)GetImageData(im_in) + MAX(0,ypos)*xdim_in 
											+ MAX(0,xpos);
	xs = MAX(0,xpos) - xpos;
	ys = MAX(0,ypos) - ypos;
	xe = MIN(xpos+xdim_out,xdim_in)-xpos;
	ye = MIN(ypos+ydim_out,ydim_in)-ypos;
	xdim_off1 = xdim_in - (xe-xs);
	xdim_off2 = xdim_out - (xe-xs);

	pim_out = (Float *)GetImageData(im_out) + ys*xdim_out + xs;

	for(y=ys; y<ye; y++,pim_in+=xdim_off1,pim_out+=xdim_off2)
	  for(x=xs; x<xe; x++,pim_in++,pim_out++)
	    *pim_out = *pim_in;
	break;
      }
      
    default:
      ERROR("Image Type not supported");
      return;
    }
  
  return;
}

/***********************************************************CommentBegin******
 *
 * -- PutSubImage -- Writes a subimage into an image
 *
 * Author :		
 *	Noel Brady
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Extracts a subimage from an image.
 * 
 * Arguments in : 	
 *	Image *im_in	: the sub-image
 *	Int xpos	: left column of the output image
 *	Int ypos	: top row of the output image
 *
 * Arguments in/out :	
 *	Image *im_out	: output image
 *
 * Arguments out :	
 *
 * Return values :	
 *	none
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	A subimage is written into the output image starting from xpos,ypos in 
 *  the output image.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
PutSubImage(Image *im_out, Image *im_in, Int xpos, Int ypos)
{
  Char *function_name = "PutSubImage";
  
  Int x,y,xdim_in,ydim_in,xdim_out,ydim_out,xdim_off;

  /*****
   *
   *	Check image dimensions and type (if im_in is >= im_out)
   *
   *****/
  if(COMPARE_IMAGE_SIZE_GOE(im_out,im_in))
    {
      ERROR_IMAGE_SIZE_GOE(im_out,im_in);
      return;
    }

  if(COMPARE_IMAGE_TYPE_EQL(im_out,im_in))
    {
      ERROR_IMAGE_TYPE_EQL(im_out,im_in);
      return;
    }

  /*****
   *
   *	Get image dimensions
   *
   *****/
  xdim_in  = GetImageSizeX(im_in);
  ydim_in  = GetImageSizeY(im_in);

  xdim_out = GetImageSizeX(im_out);
  ydim_out = GetImageSizeY(im_out);
  
  xdim_off = xdim_out - xdim_in;

  /*****
   *
   *	Work on correct image type
   *
   *****/
  switch(GetImageType(im_out))
    {
    case UCHAR_TYPE:
      {
	UChar *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *
	 *****/
	pim_out  = (UChar *)GetImageData(im_out) + ypos*xdim_out + xpos;
	pim_in = (UChar *)GetImageData(im_in);
      
	for(y=0; y<ydim_in; y++,pim_out+=xdim_off)
	  for(x=0; x<xdim_in; x++,pim_out++,pim_in++)
	    *pim_out = *pim_in;
      
	break;
      }

    case SHORT_TYPE:
      {
	SInt *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *
	 *****/
	pim_out  = (SInt *)GetImageData(im_out) + ypos*xdim_out + xpos;
	pim_in = (SInt *)GetImageData(im_in);
      
	for(y=0; y<ydim_in; y++,pim_out+=xdim_off)
	  for(x=0; x<xdim_in; x++,pim_out++,pim_in++)
	    *pim_out = *pim_in;
      
	break;
      }

    case FLOAT_TYPE:
      {
	Float *pim_in,*pim_out;
      
	/*****
	 *
	 *	    Scan image and extract subimage
	 *
	 *****/
	pim_out  = (Float *)GetImageData(im_out) + ypos*xdim_out + xpos;
	pim_in = (Float *)GetImageData(im_in);
      
	for(y=0; y<ydim_in; y++,pim_out+=xdim_off)
	  for(x=0; x<xdim_in; x++,pim_out++,pim_in++)
	    *pim_out = *pim_in;
      
	break;
      }
      
    default:
      ERROR("Image Type not supported");
      return;
    }
  
  return;
}


/***********************************************************CommentBegin******
 *
 * -- ThreshGeImage -- Binarization of an image
 *
 * Author : Michael Wollborn (TUH)		
 *	
 *
 * Created :		
 *	
 *
 * Purpose :		
 *	Binarization of an image.
 * 
 * Arguments in : 	
 *	Image *im_in	: input image
 *	Float thresh	: threshold
 *
 * Arguments in/out :	
 *	none
 *
 * Arguments out :	
 *	Image *im_out	: output image
 *
 * Return values :	
 *	none
 *
 * Side effects :	
 *	none
 *
 * Description :	
 *	The image is binarized with the threshold "thresh".
 *	All values >= thresh are set to one, all other are
 *	set to zero.
 *
 * See also :
 *	
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Void
ThreshGeImage(Image *im_in, Image *im_out, Float thresh)
{
  Char *function_name = "ThreshGeImage";
  
  Int i;

  /*****
   *
   *	Check image dimensions and type
   *
   *****/
  if(COMPARE_IMAGE_SIZE_EQL(im_in,im_out))
    {
      ERROR_IMAGE_SIZE_EQL(im_in,im_out);
      return;
    }

  if(COMPARE_IMAGE_TYPE_EQL(im_in,im_out))
    {
      ERROR_IMAGE_TYPE_EQL(im_in,im_out);
      return;
    }

  /*****
   *
   *	Work on correct image type
   *
   *****/
  switch(GetImageType(im_in))
    {
    case UCHAR_TYPE:
      {
	UChar *pim_in,*pim_out,lthresh;
	lthresh = (UChar)ROUND(thresh);
      
	/*****
	 *
	 *	    Scan image and set all values >= lthresh to "1"
	 *
	 *****/
	pim_in  = (UChar *)GetImageData(im_in);
	pim_out = (UChar *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  if(*pim_in >= lthresh)
	    *pim_out = 1;
	  else
	    *pim_out = 0;

	break;
      }

    case SHORT_TYPE:
      {
	SInt *pim_in,*pim_out,lthresh;
	lthresh = (SInt)ROUND(thresh);
      
	/*****
	 *
	 *	    Scan image and set all values >= lthresh to "1"
	 *
	 *****/
	pim_in  = (SInt *)GetImageData(im_in);
	pim_out = (SInt *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  if(*pim_in >= lthresh)
	    *pim_out = 1;
	  else
	    *pim_out = 0;

	break;
      }

    case FLOAT_TYPE:
      {
	Float *pim_in,*pim_out,lthresh;
	lthresh = thresh;
      
	/*****
	 *
	 *	    Scan image and set all values >= lthresh to "1"
	 *
	 *****/
	pim_in  = (Float *)GetImageData(im_in);
	pim_out = (Float *)GetImageData(im_out);
  
	for(i=GetImageSize(im_in); i; i--,pim_in++,pim_out++)
	  if(*pim_in >= lthresh)
	    *pim_out = 1;
	  else
	    *pim_out = 0;

	break;
      }
      
    default:
      ERROR("Image Type not supported");
      return;
    }
  
  return;
}


/*****************************************************************************/
/*********************************************** CommentBegin ****************/
/*
 *--------------------------------------------------------------
 *
 * -- CutPaste -- xxxxxxxxxxxxxx
 *
 *      Int     CutPasteImage(Image *ImaFrom,  Image *ImaTo)
 *
 * Purpose :			Get a rectangle from an image & paste into  another
 *	
 *	
 *	
 * Arguments in:
 *
 *	Image *ImaFrom
 *  long  tlx, tly		Top-left coordinates of rectangle to be copied
 *  long  brx, bry		Bottom-right coordinates of rectangle to be copied 
 *  long  x, y			Top-left coordinates where to copy the rectangle
 *
 * Arguments in/out:
 *
 *
 *
 * Arguments out:
 *
 *	Image *ImaTo
 *
 *
 * Description:			Works on SHORT_TYPE, FLOAT_TYPE and UCHAR_TYPE images
 *
 *
 * Return values:		1 on error, otherwise 0 
 *
 *
 * Side effects:		None
 *
 *
 * See also:
 *
 *
 *--------------------------------------------------------------
 */
/*************************************************CommentEnd*******************/
Int CutPasteImage (Image *ImaFrom, Int tlx, Int tly,
			                  Int brx, Int bry,
			  Image *ImaTo,   Int x,   Int y)
{
	 SInt  	*pIns,  *pOuts;
	 Float	*pInf,  *pOutf;
	 UChar 	*pInuc, *pOutuc;
	 Int	dim_h, dim_v;		/* Must be signed! */
	 UInt	i, j;

	 if (tlx < 0 || brx < 0 || tly < 0 || bry < 0)
	 {
		  runerr("ERROR(cutpaste): invalid input rectangle sizes.", -1);
	 }

	 dim_h = brx - tlx;
	 dim_v = bry - tly;

	 if (dim_h < 0 || dim_v < 0 || dim_h >= GetImageSizeX(ImaFrom) ||
		                           dim_v >= GetImageSizeY(ImaFrom))
	 {
		  runerr("ERROR(cutpaste): invalid input rectangle sizes.", -1);
	 }
	 
	 if (x < 0 || y < 0 || x + dim_h >= GetImageSizeX(ImaTo) ||
		                   y + dim_v >= GetImageSizeY(ImaTo))
	 {
		  runerr("ERROR(cutpaste): invalid output rectangle coords.", -1);
	 }

	 switch(GetImageType(ImaFrom))
	 {
	   case SHORT_TYPE:
		  if (GetImageType(ImaTo) != SHORT_TYPE)
		  {
			   runerr("ERROR(CutPaste): image not of expected type",1);
		  }

		  /* Copy the data from Imafrom in ImaTo line by line, */
		  /* starting from the top							  */
		  /* ------------------------------------------------- */

		  for (i = 0; i <= dim_v; i++)
		  {
			   /* Get the address of the left-most first pixel */
			   /* of the line to be copied                     */
			   /* -------------------------------------------- */
			   pIns = (SInt *) GetImageData(ImaFrom) +
			                   GetImageSizeX(ImaFrom)*(tly+i)+tlx;

			   /* Set pointer in the output image */
			   /* ------------------------------- */
			   pOuts = (SInt *) GetImageData(ImaTo) + 
			                    GetImageSizeX(ImaTo)*(y+i) + x;

			   /* Copy the current line */
			   /* --------------------- */
			   for (j = 0; j <= dim_h; j++)
			   {
					*pOuts++ = *pIns++;
			   }
		  }
		  break;

	   case FLOAT_TYPE:
		  if (GetImageType(ImaTo) != FLOAT_TYPE)
		  {
			   runerr("ERROR(CutPaste): image not of expected type",1);
		  }

		  /* Copy the data from Imafrom in ImaTo line by line, */
		  /* starting from the top							  */
		  /* ------------------------------------------------- */

		  for (i = 0; i <= dim_v; i++)
		  {
			   /* Get the address of the left-most first pixel */
			   /* of the line to be copied                     */
			   /* -------------------------------------------- */
			   pInf  = (Float *) GetImageData(ImaFrom) + 
			                     GetImageSizeX(ImaFrom)*(tly+i)+tlx;

			   /* Set pointer in the output image */
			   /* ------------------------------- */
			   pOutf = (Float *) GetImageData(ImaTo) + 
			                     GetImageSizeX(ImaTo)*(y+i) + x;

			   /* Copy the current line */
			   /* --------------------- */
			   for (j = 0; j <= dim_h; j++)
			   {
					*pOutf++ = *pInf++;
			   }
		  }
		  break;

	   case UCHAR_TYPE:
		  if (GetImageType(ImaTo) != UCHAR_TYPE)
		  {
			   runerr("ERROR(CutPaste): image not of expected type",1);
		  }

		  /* Copy the data from Imafrom in ImaTo line by line, */
		  /* starting from the top							  */
		  /* ------------------------------------------------- */

		  for (i = 0; i <= dim_v; i++)
		  {
			   /* Get the address of the left-most first pixel */
			   /* of the line to be copied                     */
			   /* -------------------------------------------- */
			   pInuc  = (UChar *) GetImageData(ImaFrom) + 
			                      GetImageSizeX(ImaFrom)*(tly+i)+tlx;

			   /* Set pointer in the output image */
			   /* ------------------------------- */
			   pOutuc = (UChar *) GetImageData(ImaTo) + 
			                      GetImageSizeX(ImaTo)*(y+i) + x;

			   /* Copy the current line */
			   /* --------------------- */
			   for (j = 0; j <= dim_h; j++)
			   {
					*pOutuc++ = *pInuc++;
			   }
		  }
		  break;

	   default:
		  printf("Image type >>%d<< not supported\n",GetImageType(ImaFrom));
		  runerr("ERROR(CutPaste) ::: image not of expected type",1);
	 }

	 return 0;
}  	/* CutPasteImage */


    
/***********************************************************CommentBegin******
 *
 * -- ValidCoordinate -- 
 *
 * Author :		
 *	Noel O'Connor
 *
 * Created :		
 *	3-March-1996 
 *
 * Purpose :		
 *	This function checks if a pixel coordinate is within an image
 * 
 * Arguments in : 	
 *	Int x - x-coordinate
 *	Int y - y-coordinate
 *	Int dim_x - horizontal dimension of image
 *	Int dim_y - vertical dimension of image
 *
 * Arguments in/out :	
 *	-
 *
 * Arguments out :	
 *	-
 *
 * Return values :	
 *	TRUE if coordinate is in image, FALSE otherwise
 *
 * Side effects :	
 *	-
 *
 * Description :	
 *	-
 *
 * See also :
 *	-
 *
 * Modified :		
 *	
 *
 ***********************************************************CommentEnd********/

Int
ValidCoordinate(Int x, Int y, Int dim_x, Int dim_y)
{
  if((x>=0) && (y>=0) && (x<dim_x) && (y<dim_y))
    return(TRUE);
  else
    return(FALSE);
}


/***********************************************************CommentBegin******
 *
 * -- MaskImage -- 
 *
 *			Int MaskImage(Image *im_in, Image *mask1, Image *mask2, Image *im_out)
 *
 * Author :		
 *	Noel Brady
 *
 * Created :		
 *	6-6-96 
 *
 * Purpose :		
 *	Mask out the pixels in an image
 * 
 * Arguments in : 	
 *	Image *im_in - input image
 *	Image *mask - mask image
 *
 * Arguments in/out :	
 *	-
 *
 * Arguments out :	
 *	Image *im_out - output masked image
 *
 * Return values :	
 *	The number of pixels in the mask image != 0.
 *
 * Side effects :	
 *	-
 *
 * Description :	
 *	-
 *
 * See also :
 *	-
 *
 * Modified :		
 *	2-Dec-96 Jan De Lameillieure (HHI) : calculate SNR only in pels that belong to 
 *                                           original AND reconstructed segment mask
 *
 ***********************************************************CommentEnd********/

Int MaskImage(Image *im_in, Image *mask1, Image *mask2, Image *im_out)
{
	switch(GetImageType(im_in))
    {
    case SHORT_TYPE:
      return MaskImageI(im_in,mask1,mask2,im_out);
      break;

    default:
      printf("MaskImage: Image type >>%d<< not supported\n", GetImageType(im_in));
      return -1;
    }
}


Int MaskImageI(Image *im_in, Image *mask1, Image *mask2, Image *im_out)
{
	SInt *pi, *po, *pm1, *pm2, *pe;

	Int n;

  Char *function_name = "MaskImage";

	switch(GetImageType(im_in))
    {
    case SHORT_TYPE:
      break;

    default:
      printf("MaskImage: Image type >>%d<< not supported\n", GetImageType(im_in));
      return -1;
    }

	switch(GetImageType(mask1))
    {
    case SHORT_TYPE:
      break;

    default:
      printf("MaskImage: Image type >>%d<< not supported\n", GetImageType(mask1));
      return -1;
    }

	switch(GetImageType(mask2))
    {
    case SHORT_TYPE:
      break;

    default:
      printf("MaskImage: Image type >>%d<< not supported\n", GetImageType(mask2));
      return -1;
    }

	switch(GetImageType(im_out))
    {
    case SHORT_TYPE:
      break;

    default:
      printf("MaskImage: Image type >>%d<< not supported\n", GetImageType(im_out));
      return -1;
    }

  if(COMPARE_IMAGE_SIZE_EQL(im_in,mask1))
    {
      ERROR_IMAGE_SIZE_EQL(im_in,mask1);
      return -1;
    }

  if(COMPARE_IMAGE_SIZE_EQL(im_in,mask2))
    {
      ERROR_IMAGE_SIZE_EQL(im_in,mask2);
      return -1;
    }

  if(COMPARE_IMAGE_SIZE_EQL(mask1,im_out))
    {
      ERROR_IMAGE_SIZE_EQL(mask1,im_out);
      return -1;
    }

	pi = (SInt *) GetImageData(im_in);
	po = (SInt *) GetImageData(im_out);
	pm1 = (SInt *) GetImageData(mask1);
	pm2 = (SInt *) GetImageData(mask2);
	pe = po + GetImageSize(im_out);

	SetConstantImage(im_out,0);
	n = 0;
	while (po != pe) {
		if ((*pm1) && (*pm2))
			{
				*po = *pi;n++;
			}
		pi++;po++;pm1++;pm2++;
	}

	return n;
}


