/*----------------------- Rogaland University Centre -------------------------
## \section{qmf (Quadrature Mirror Filterbank)}
## 
## \subsection*{Purpose}
## Filter bank based on tree-structured prototype filters.
## Employing both decimated and full-rate analysis filterbank
## and decimated synthesis filterbank, and both general FIR
## filters and IIR filters based on polyphase allpass filters.
## 
## \subsection*{Description}
## The program will read the filterparameters from files specified
## in the parameters and perform the filtering on the input image.
## 
## Note that a filter identifier starting with `f' is assumed to
## represent a FIR filter, an identifier starting with `i'
## is a polyphase-allpass IIR filter and an identifier starting
## with `w' is a wavelet FIR filter.
## 
## \subsection*{Usage }
## \begin{verbatim}
## qmf [-sn struct<integer> | 
##      -s struct<string> | 
##      -sb struct<string>]
##     [-f filter-id.<string>] [-d] 
##     [-analysis | -synthesis]
##     [imagefile]
## \end{verbatim}
## \begin{description}
##   \item[-f]
## 	Filter specification file to be used as default filter [f16b]
##   \item[-s]
## 	Filter structure string, defined as a quad tree.
## 	See the module {\em filtertree} (CompileFiltStructTree)
## 	for a description of the structure.
##   \item[-sb]
## 	Filterstructure string, defined as a binary tree.
## 	See the module {\em filtertree} (CompileFiltStructTreeFromBin)
## 	for a description of the structure.
##   \item[-sn]
## 	Filter structure number, between 1 and 11 (default 4).
## 	See the module {\em filtertree} (BuildFiltStructTree)
## 	for a description of the structure.
##   \item[-d]
## 	Critically decimated filter (default FALSE).
## 	When a decimated filter is used, there will be one output frame
## 	per input frame (with the same size as the input --- a mosaic
## 	of the decimated subbands).
## 	When a fullrate filter is used, there will be one output
## 	frame per subband per input frame --- the same size as the 
## 	input image.
##   \item[-analysis]
## 	Analysis filterbank (default TRUE).
##   \item[-synthesis]
## 	Synthesis filterbank.
## \end{description}
## 
## 
## 
## 
## 
## \subsection*{Date created}
## March 1993
 *------------------------ Rogaland University Centre -----------------------*/



/* Header files */
#include "qmf.h"

static Flag_Format flagfmt[]={
  {"sn",{"s","sb",LASTFLAG},1,{
    {PTINT,"4","struct"},
    LASTPARAMETER}},
  {"f",{LASTFLAG},1,{
    {PTSTRING,"f16b","filter-id."},
    LASTPARAMETER}},
  {"sb",{"s","sn",LASTFLAG},1,{
    {PTBOOLEAN,"FALSE"},
    {PTSTRING,"(,)","struct"},
    LASTPARAMETER}},
  {"s",{"sb","sn",LASTFLAG},1,{
    {PTBOOLEAN,"FALSE"},
    {PTSTRING,"(,)","struct"},
    LASTPARAMETER}},
  {"d",{LASTFLAG},0,{
    {PTBOOLEAN,"FALSE"}, 
    LASTPARAMETER}},
  {"analysis",{"synthesis",LASTFLAG},0,{
    {PTBOOLEAN,"TRUE"},
    LASTPARAMETER}},
  {"synthesis",{"analysis","d",LASTFLAG},0,{
    {PTBOOLEAN,"FALSE"},
    LASTPARAMETER}},
  LASTFLAG};

int types[]={PFFLOAT,LASTTYPE};

main(argc,argv)
     int argc;
     char *argv[];
{
  FILE *fp;
  Filename filename;		/* Input image file */
  Filename filter;		/* Filterfile name */
  struct header hd;		/* Input file header */
  struct header hdp;		/* Converted input header */
  struct header hd_c;		/* Byte-formatted output header */
  int ro,kappa,i,j,m,n,k,fr;
  int argcc;
  int method, method_c;		/* Conversion method */
  float **y;			/* Pointer to the image pixels */
  int  struct_no;		/* Subband split structure number */
  Boolean decimateFlag;		/* Decimate the output */
  Boolean structStrFlag;	/* Structure string specified (quad) */
  Boolean structStr1DFlag;	/* Structure string specified (bin) */
  Boolean analysisFlag;		/* Analysis filtering */
  Boolean synthesisFlag;	/* Synthesis filtering */
  char *structStr;		/* Structure description string (quad) */
  char *structStr1D;		/* Structure description string (bin) */
  FilterTree filtertree;	/* Filter structure tree */
  int n_subbnds;		/* Number of subbands */

  /* Initialize variables */
  Progname=strsave(*argv);
  parseargs(argc,argv,flagfmt,&struct_no,&filter,
	    &structStr1DFlag,&structStr1D,&structStrFlag,&structStr,
	    &decimateFlag,&analysisFlag,&synthesisFlag,
	    FFONE,&filename);
 
  /* Read the input image and check the dimensions */
  fp = hfopenr(filename);
  fread_hdr_a(fp,&hd,filename);
  method = fset_conversion(&hd,&hdp,types,filename);

  /* Build the filter structure binary tree */
  if (structStrFlag) 
    filtertree = CompileFiltStructTree(structStr,filter);
  else if (structStr1DFlag) 
    filtertree = CompileFiltStructTreeFromBin(structStr1D,filter);
  else
    filtertree = BuildFiltStructTree(struct_no,filter);
  n_subbnds = getNSubbands (filtertree);
  PrintFiltStructTree (filtertree);
  
  /* Check the image dimension against the filter structure */
  if ( ((hdp.rows % ((int) pow(2.0,getNSubLevels(filtertree)+1.0))) !=0) 
      || ((hdp.cols % ((int) pow(2.0,getNSubLevels(filtertree)+1.0))) !=0)) {
    fprintf(stderr,"The image dimensions are not compatible with");
    fprintf(stderr," the filterstructure\n");
    exit(1);
  }
  
  /* Write the output image header */
  if (synthesisFlag) {
    /* Convert the output to byte format */
    method_c=find_method(PFBYTE,PFFLOAT);
    dup_header(&hd,&hd_c);
    setformat(&hd_c,PFBYTE);
    write_headeru(&hd_c,argc,argv);
  }
  else {
    if (!decimateFlag)
      /* If the output is not to be decimated, then each subband will
	 be in a separate frame */
      hdp.num_frame *= n_subbnds;

    /* Write the output as a floating point image */
    write_headeru(&hdp,argc,argv);
  }
  
  /* Perform the filtering on each input frame */
  for (fr=0;fr<hd.num_frame;fr++) {
    /* Read the input frame */
    fread_imagec(fp,&hd,&hdp,method,fr,filename);
    y=(float **) matrix2(&hdp);
    
    /* Filter the image, and write the result */
    if (synthesisFlag) {	/* Synthesis filter bank */
      dec_synthesis_bank(&hdp,y,filtertree,fr);
    }
    else {	/* Analysis filter bank */
      h_linscale(&hdp,&hdp,1.0,-128.0);
      if (decimateFlag)
	dec_analysis_bank(&hdp,y,filtertree,fr);
      else
	analysis_bank (y,filtertree,hdp.rows,hdp.cols);
    }
  
    fprintf(stderr,"\n");
    free_matrix2(y);
  }
}