factor.c

FIR2LIFT is a program to factor a wavelet FIR filter pair into lifting steps.

/*
FILE : FACTOR.C

Functions to factor a polyphase matrix into lifting steps.

(C) C. Valens

Created     : 08/09/1999
Last update : 24/09/1999
*/


#include "redir.h"
#include "zpoly.h"
#include "polyphas.h"
#include "gcd.h"
#include "factor.h"

#include <stdlib.h>


/*
 * Local prototypes.
 */
factor_list *factor_list_goto_end(factor_list *fl);
factor_list *factor_list_remove_multiples(factor_list *fl);
void factor_list_write_one(factor_list *fl);
void factor_list_write_list(factor_list *fl);


/*
 * Debug varable.
 */
static int factor_list_count = 0;

/*
 * Counter for factorizations.
 */
int factorizations=1;

factor_list *factor_list_create(void)
{
  factor_list *l;
  if ((l=malloc(sizeof(factor_list)))==NULL) {
    redir_printf("factor_list_create() : ");
    zpoly_error(1);
  }
  memset(l,0,sizeof(factor_list));

  l->tag = factor_list_count;
  factor_list_count++;

  return l;
}


void factor_list_destroy_one(factor_list *fl)
{
  factor_list *l1;

  l1 = fl;

  if (l1) {
    polyphase_matrix_destroy(l1->P);
    polyphase_matrix_destroy(l1->F);
    free(l1);
  }
}


void factor_list_destroy(factor_list *fl)
{
  factor_list *l1, *l2;

  l1 = fl;

  while (l1) {
    polyphase_matrix_destroy(l1->P);
    polyphase_matrix_destroy(l1->F);
    l2 = l1;
    l1 = l1->next;
    free(l2);
  }
}


/*
 * Creates a copy of a factor list element.
 */
factor_list *factor_list_copy(factor_list *fl)
{
  factor_list *l = fl;
  factor_list *temp;
  factor_list *flnew = NULL;

  while (l) {

    temp = factor_list_create();
    if (flnew==NULL) {
      flnew = temp;
    }
    temp->P = polyphase_matrix_copy(fl->P);
    temp->F = polyphase_matrix_copy(fl->F);
    temp = temp->next;

    l = l->next;
  }
  return flnew;
}


/*
 * Returns a pointer to the last element of a factor list.
 */
factor_list *factor_list_goto_end(factor_list *fl)
{
  factor_list *l, *last = NULL;
  l = fl;
  while (l) {
    last = l;
    l = l->next;
  }
  return last;
}


/*
 * Remove all multiple entries from a factor list.
 */
factor_list *factor_list_remove_multiples(factor_list *fl)
{
  factor_list *l1, *l2, *l3;

  l1 = fl;

//  redir_printf("{multiple list}\n");
//  factor_list_write_list(l1);
//  redir_printf("{multiple list end}\n");

  while (l1) {
    l2 = l1->next;
    l3 = l1;
    while (l2) {

//      factor_list_write_one(l2);

      if ((polyphase_matrix_equal(l1->P,l2->P)) &&
          (polyphase_matrix_equal(l1->F,l2->F))) {

//        redir_printf("{l1==l2}\n");
//        factor_list_write_one(l1);
//        factor_list_write_one(l2);
//        factor_list_write_one(l3);

        l3->next = l2->next;
        factor_list_destroy_one(l2);
        l2 = l3;

//        factor_list_write_one(l3->next);
//        redir_printf("{end}\n");
      }
      l3 = l2;
      l2 = l2->next;
    }
    l1 = l1->next;
  }

  return fl;
}


/*
 * Debug function.
 *
 * Writes an element of a factor tree/list to the output stream.
 */
void factor_list_write_one(factor_list *fl)
{
  if (fl) {
    polyphase_matrix_write(fl->P);
    polyphase_matrix_write(fl->F);
  }
}


/*
 * Debug function.
 *
 * Writes a list of factors to the output stream.
 */
void factor_list_write_list(factor_list *fl)
{
  factor_list *l = fl;

  while (l) {
    factor_list_write_one(l);
    l = l->next;
  }
}


/*
 * Debug function.
 *
 * Writes all elements in a factor tree to the output stream, not in
 * any meaningfull order.
 */
void factor_list_write_tree(factor_list *fl)
{
  factor_list *l = fl;

  if (l) {
    factor_list_write_tree(l->lift);
    factor_list_write_tree(l->next);
  }
  factor_list_write_one(l);
}


/*
 * Internal function.
 *
 * Assemble a lifting step.
 */
factor_list *factor_build(factor_list *fl, zpoly *p1, zpoly *p2,
                          zpoly *p3, zpoly *p4, zpoly *p5, zpoly *p6)
{
//  redir_printf("~~~\n");
  fl->P.row1.p1 = zpoly_copy(p1);
  fl->P.row1.p2 = zpoly_copy(p2);
  fl->P.row2.p1 = zpoly_copy(p3);
  fl->P.row2.p2 = zpoly_copy(p4);
  fl->F.row1.p1 = zpoly_create(0,1.0);
  if (p5) fl->F.row1.p2 = zpoly_copy(p5);
  else fl->F.row1.p2 = NULL;
  if (p6) fl->F.row2.p1 = zpoly_copy(p6);
  else fl->F.row2.p1 = NULL;
  fl->F.row2.p2 = zpoly_create(0,1.0);

//  factor_list_write_one(fl);
//  redir_printf("+++\n");
  return fl;
}


/*
 * Extract a primal lifting step (the upper triangular matrix):
 *
 * /        \    /      \ /      \
 * | he  ho | => | a  b | | 1  L | => he=a, ho=a+bL
 * | ge  go |    | c  d | | 0  1 |    ge=c, go=c+dL
 * \        /    \      / \      /
 *
 * PRIMAL1 starts with (he,ho), PRIMAL2 starts with (ge,go)
 *
 * OR
 *
 * Extract a dual lifting step (the lower triangular matrix):
 *
 * /        \    /      \ /      \
 * | he  ho | => | a  b | | 1  0 | => he=a+bL, ho=b
 * | ge  go |    | c  d | | L  1 |    ge=a+cL, go=d
 * \        /    \      / \      /
 *
 * DUAL1 starts with (he,ho), DUAL2 starts with (ge,go)
 *
 * Returns a list with all possible values for L.
 */
factor_list *factor_extract_step(zpoly *p1, zpoly *p2, zpoly *p3, zpoly *p4,
                                 int type)
{
  zpoly *p, *temp, *div, *mod;
  gcd_list *gcds, *gcd;
  factor_list *r, *result = NULL;

  gcds = gcd_euclides(p1,p2);

  gcd = gcds;
  while (gcd) {

    div = gcd->divmod.p1;
    mod = gcd->divmod.p2;

    /*
     * Avoid trivial solutions and make sure that degrees diminish.
     */
    if ((!zpoly_equals_zero(div)) &&
        (zpoly_degree(mod)<=zpoly_degree(p1))) {

      p = zpoly_mul(p4,div);
      temp = zpoly_sub(p3,p);
      if (zpoly_degree(temp)<=zpoly_degree(p3)) {

        if (result==NULL) {
          result = factor_list_create();
          r = result;
        }
        else {
          r->next = factor_list_create();
          r = r->next;
        }

        /*
         * Select the correct way to assemble the matrices.
         */
        switch (type) {

          case PRIMAL1:
            factor_build(r,p2,mod,p4,temp,div,NULL);
          break;

          case PRIMAL2:
            factor_build(r,p4,temp,p2,mod,div,NULL);
          break;

          case DUAL1:
            factor_build(r,mod,p2,temp,p4,NULL,div);
          break;

          case DUAL2:
            factor_build(r,temp,p4,mod,p2,NULL,div);
          break;
        }
      }
      zpoly_destroy(p);
      zpoly_destroy(temp);
    }

    gcd = gcd->next;
  }

  gcd_destroy_list(gcds);

  return result;
}



factor_list *factor_extract_lifting_step(polyphase_matrix m, int type)
{
  factor_list *r, *last, *result = NULL;

  if (type==PRIMAL) {
    result = factor_extract_step(m.row1.p2,m.row1.p1,m.row2.p2,m.row2.p1,
                                 PRIMAL1);

    last = factor_list_goto_end(result);

    r = factor_extract_step(m.row2.p2,m.row2.p1,m.row1.p2,m.row1.p1,
                            PRIMAL2);

    if (result==NULL) result = r;
    else last->next = r;
  }
  else {
    result = factor_extract_step(m.row1.p1,m.row1.p2,m.row2.p1,m.row2.p2,
                                 DUAL1);

    last = factor_list_goto_end(result);

    r = factor_extract_step(m.row2.p1,m.row2.p2,m.row1.p1,m.row1.p2,
                            DUAL2);

    if (result==NULL) result = r;
    else last->next = r;
  }

  result = factor_list_remove_multiples(result);

  return result;
}


/*
 * Expands a matrix into lifting steps, ie alternating dual & primal
 * matrices.
 */
factor_list *factor_extract_lifting_steps(factor_list *fl,int type)
{
  factor_list *l = fl;

  while (l) {

    /*
     * Do a primal or dual lifting step.
     */
    l->lift = factor_extract_lifting_step(l->P,type);

    if (type==PRIMAL) {
      /*
       * Continue with a dual lifting step.
       */
      factor_extract_lifting_steps(l->lift,DUAL);
    }
    else {
      /*
       * Continue with a primal lifting step.
       */
      factor_extract_lifting_steps(l->lift,PRIMAL);
    }
    l = l->next;
  }

  return fl;
}


/*
 * Stack for internal use only.
 */
static factor_list *stack_pointer = NULL;

void factor_list_push(factor_list *fl)
{
  factor_list *temp = factor_list_copy(fl);

  if (stack_pointer==NULL) {
    stack_pointer = temp;
  }
  else {
    stack_pointer->next = temp;
    temp->prev = stack_pointer;
    stack_pointer = stack_pointer->next;
  }
}

void factor_list_pop(void)
{
  factor_list *temp;

  if (stack_pointer) {

    temp = stack_pointer->prev;
    factor_list_destroy(stack_pointer);
    stack_pointer = temp;
    if (stack_pointer) {
      stack_pointer->next = NULL;
    }
  }
}


factor_list *factor_list_dump_stack(void)
{
  factor_list *temp = stack_pointer;
  factor_list *result = NULL;
  factor_list *l;

  while (temp) {

    if (result==NULL) {
      result = factor_list_create();
      l = result;
    }
    else {
      l->lift = factor_list_create();
      l = l->lift;
    }	
    l->P = polyphase_matrix_copy(temp->P);
    l->F = polyphase_matrix_copy(temp->F);

    temp = temp->prev;
  }
  
  return result;
}


/*
 * Simplify life by globally fixing the list.
 */
static factor_list *factor_root = NULL;
static factor_list *last = NULL;

/*
 * Returns a list of all possible factorizations.
 */
factor_list *factor_collect(factor_list *fl)
{
  if (fl) {

    factor_list_push(fl);

    factor_collect(fl->lift);

    /*
     * Build a double linked list of all possible factorizations.
     */
    if (factor_root==NULL) {
      factor_root = factor_list_dump_stack();
      last = factor_root;
    }
    else {
      last->next = factor_list_dump_stack();
      if (last->next) last->next->prev = last;
      last = last->next;
    }

    factor_list_pop();

    if (fl->next) {
      factor_collect(fl->next);
    }
  }
  return factor_root;
}


/*
 * Takes a (wavelet) filter pair as input and returns a list of all
 * possible factorizations in lifting steps.
 */
factor_list *factor_this_baby_now(zpoly_pair filters)
{
  factor_list *fl1, *fl2, *result=NULL;

  fl1 = factor_list_create();
  fl1->P = polyphase_matrix_build(filters);
  fl1->F.row1.p1 = zpoly_create(0,1);
  fl1->F.row2.p2 = zpoly_create(0,1);
  fl2 = factor_list_copy(fl1);

//  factor_list_write_one(fl1);

  fl1 = factor_extract_lifting_steps(fl1,PRIMAL);

//  factor_list_write_one(fl2);

  fl2 = factor_extract_lifting_steps(fl2,DUAL);
  factor_collect(fl1);
  result = factor_collect(fl2);

/*
 * Hier moet nog het een en ander komen, denk ik.
 */

  factor_list_destroy(fl1);
  factor_list_destroy(fl2);

  return result;
}


void gimme_my_lifting_steps(factor_list *fl)
{
  factor_list *l1 = fl;
  factor_list *l2;
  zpoly *d;

  if (l1) {
    d = polyphase_matrix_determinant(l1->P);
    redir_printf("*** determinant = [");
    zpoly_write(d);
    redir_printf("]\n\n");
  }

  while (l1) {
    l2 = l1;

    if ((l2) && (l1->next)) {
      redir_printf("*** %d:\n",factorizations);
      polyphase_matrix_write(l2->P);
      factorizations++;
    }
    while (l2) {
      if (l2->lift) {
        polyphase_matrix_write(l2->F);
      }
      l2 = l2->lift;
    }

    if (l1->next) redir_printf("\n");

    l1 = l1->next;
  }

  factorizations--;
  zpoly_destroy(d);
}