/* 
   Copyright (C) 1999 E. H. Haley

   wtd does the 2d 3color wavelet transform, but in each row at a time
   all the way down then each column.  This produces quantitatively
   different results than wtliketargd which emulates the one level at
   a time procedure of targa.c (was targd.c), and which should be
   used.

 */

#include <stdlib.h>
#include "wto.h"
#include "daub4.h"

void wtd(float a[], int w, int h, int c,
	 int isign, void (*rowtxfn)(float [], int, int))
{ 
  int i,j,k,nt;  float *wksp;
  
  wksp = (float *) malloc ((w>h?w:h)*sizeof(float));
  for (k=0; k<c; k++)
    {
      if (w>=4) //do rows
	for (j=0; j<h; j++)
	  {
	    for (i=0; i<w; i++)
	      wksp[i]=a[c*(w*j+i)+k];
	    if (isign >= 0)
	      for(nt=w;nt>=4;nt >>= 1)
		//		(*rowtxfn)(wksp-1,nt,isign);
		(*rowtxfn)(wksp,nt,isign);
	    else
	      for(nt=4;nt<=w;nt <<= 1)
		//		(*rowtxfn)(wksp-1,nt,isign);
		(*rowtxfn)(wksp,nt,isign);
	    for (i=0; i<w; i++)
	      a[c*(w*j+i)+k]=wksp[i];
	  }
      
      if(h>=4) // do cols
	for (i=0; i<w; i++)
	  {
	    for (j=0; j<h; j++)
	      wksp[j]=a[c*(w*j+i)+k];
	    if (isign >= 0)
	      for(nt=h;nt>=4;nt >>= 1)
		(*rowtxfn)(wksp-1,nt,isign);
	    else
	      for(nt=4;nt<=h;nt <<= 1)
		(*rowtxfn)(wksp-1,nt,isign);
	    for (j=0; j<h; j++)
	      a[c*(w*j+i)+k]=wksp[j];
	  }
    }
  free(wksp);
}

void wtliketargd(float a[], int w, int h, int c,
		 int isign, void (*rowtxfn)(float [], int, int))
{
  int i,j,k,nt;  float *wksp;

  wksp = (float *) malloc ((w>h?w:h)*sizeof(float));

  for(nt= (isign>0?w:4);  nt>= (isign>0?4:w);  nt = (isign>0)?nt>>1:nt<<1 )
    //requires w==h
    {
      for(k=0; k<c; k++)
	{
	  for(j=0; j<nt; j++)
	    {
	      for(i=0; i<nt; i++)
		wksp[i]=a[c*(w*j+i)+k];
	      //	      (*rowtxfn)(wksp-1,nt,isign);
	      (*rowtxfn)(wksp,nt,isign);
	      for(i=0; i<nt; i++)
		a[c*(w*j+i)+k]=wksp[i];
	    }
	  for(i=0; i<nt; i++)
	    {
	      for(j=0; j<nt; j++)
		wksp[j]=a[c*(w*j+i)+k];
	      //	      (*rowtxfn)(wksp-1,nt,isign);
	      (*rowtxfn)(wksp,nt,isign);
	      for(j=0; j<nt; j++)
		a[c*(w*j+i)+k]=wksp[j];
	    }
	}
    }
  free(wksp);
}