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

   args jpeg file to read from, datafile to append to.

   Records (filename, size, and) top 40 DAUB4 wavelet coeffs for the
   whole target image, and each quadrant, each sixteenth, sixtyfourth,
   etc., only it does it in reverse order.  It's yanking intermediate
   arrays from the one wavelet transforming process.

   wtd() is faster _because_ it doesn't do it one layer at a time, but
   it doesn't produce the same results (off by one part in 50,
   sometimes, at least), so tirade.c uses wtliketargd(), of all
   things, which uses, well, you know.

   Note "tree" is so-called, but it's an array of _res_ not _node_ (in
   a 4-ary tree configuration, children of k at 4k+1, 4k+2, 4k+3,
   4k+4).  In snatches.c, they'll be read into a res and copied
   backwards into another treeLike array of res, with a similar array
   [2][targ.nnodes] of type tile, that holds their contenders for
   matches.

   targa extends res to include the actual average luminance, not with
   these funky factors, for each target subsquare and each tile.

   Getting these average Ys after YIQ but before the WT.  Hey, should
   just calc them once for the 16x16's, then get next largers by
   averaging their littlers.  parent=child/4?

*/

#include <stdio.h>
#include <jpeglib.h>
#include "lib/loadj.h"
#include "lib/top.h"
#include "lib/wto.h"
#include "lib/color.h"
#include "lib/tile.h"
#include "lib/daub4.h"

int main(int aardc, char **aardv)
{
  int i,j,k,l,p, nt, w,h,c, subw,subh, nw, nh, x,y, node, ind,outd;
  JSAMPLE *arr;
  float heap[40], *scr, *scr2, *line;
  unsigned long index[40];
  FILE *datfile;
  tfile curimg;
  res *tree,*treerev;

  if (aardc<3) exit(0);
  if((arr=rjf(*++aardv, &w,&h,&c))==NULL)
    { fprintf(stderr, "targa: quitting\n"); exit(0); }

  strcpy(curimg.fname,*aardv);
  curimg.side=w;

  if ((datfile=fopen(*++aardv,"w"))==NULL)
    { fprintf(stderr, "targa: can't open datfile %s\n", *aardv); exit(1); }

  scr=(float *)malloc(w*h*c* sizeof(float));
  scr2=(float *)malloc(w*h*c* sizeof(float));
  tree=(res *)malloc((w*h/32)*sizeof(res)); /* /32 not /16*3, to be safe */
  line=(float *)malloc((w>h?w:h)*sizeof(float));

  for(i=0; i<w*h; i++)
    rgbk256yiqy(arr+c*i,scr+c*i,c);

  //---------------------------------------------  

  subw=16; subh=16; nw=w/16; nh=h/16; node=0;  
  for(nt=w;nt>=16;nt >>= 1)  //requires w==h again
    {
      //nsquares=nw*nh; side=subw,subh
      for(x=0; x<nw; x++)
	for(y=0; y<nh; y++) //for each subsquare of image
	  {
	    for(i=0; i<subw; i++)
	      for(j=0; j<subh; j++)
		{
		  tree[node].avgy
		    += scr[c*(w*(j+y*subh)+(i+x*subw))+0]; //k=0
		}
	    tree[node].avgy/=(float)(subw*subh);
	    node++;
	  }
      subw*=2; subh*=2; nw/=2; nh/=2;
    }
  
  //-------------------------------------------------  
  subw=4; subh=4; nw=w/4; nh=h/4; node=0;
  
  //code for wtliketargd() was lifted straight from this, thence the name:
  for(nt=w;nt>=4;nt >>= 1)  //requires w==h again
    {
      for(k=0; k<c; k++)
	{
	  for(j=0; j<nt; j++)
	    {
	      for(i=0; i<nt; i++)
		line[i]=scr[c*(w*j+i)+k];
	      //	      daub4(line-1,nt,1);
	      dobb4(line,nt,1);
	      for(i=0; i<nt; i++)
		scr[c*(w*j+i)+k]=line[i];
	    }
	  for(i=0; i<nt; i++)
	    {
	      for(j=0; j<nt; j++)
		line[j]=scr[c*(w*j+i)+k];
	      //	      daub4(line-1,nt,1);
	      dobb4(line,nt,1);
	      for(j=0; j<nt; j++)
		scr[c*(w*j+i)+k]=line[j];
	    }
	}
      if(subw>=16) //record them
	{
	  //nsubsquares = nw*nh
	  //subsquare's sides = subw,subh
	  //subw*nw=w, subh*nh=h
	  //BTW, always: nw=nt/4
	  //
	  for(x=0; x<nw; x++)
	    for(y=0; y<nh; y++) //for each subsquare of image
	      {
		for(k=0; k<c; k++)
		  {
		    for(l=2; l<=subw/2; l<<=1)
		      for(i=0; i<subw/l; i++)
			for(j=0; j<subh/l; j++)
			  {
			    scr2[c*(subw*(j+subh/l)+(i+subw/l))+k]
			      = scr[c*(w*(j+y*subh/l+h/l)+(i+x*subw/l+w/l))+k];
			    scr2[c*(subw*(j       )+(i+subw/l))+k]
			      = scr[c*(w*(j+y*subh/l    )+(i+x*subw/l+w/l))+k];
			    scr2[c*(subw*(j+subh/l)+(i       ))+k]
			      = scr[c*(w*(j+y*subh/l+h/l)+(i+x*subw/l    ))+k];
			  }
		    for(i=0; i<2; i++)
		      for(j=0; j<2; j++)
			scr2[c*(subw*j+i)+k] 
			  = scr[c*(w*(j+y*2)+(i+x*2))+k];
		  }
		
		top(40,subw*subh*c,scr2-1,heap-1,index-1);
		for (i=0; i<40; i++)
		  {
		    tree[node].side = subw;
		    tree[node].entry[i].index[0] = (index[i]/c)%subw;
		    tree[node].entry[i].index[1] = index[i]/(c*subw);
		    tree[node].entry[i].index[2] = index[i]%c;
		    tree[node].entry[i].val = heap[i];
		  }
		node++;
	      }
	}
      subw*=2; subh*=2; nw/=2; nh/=2;
    }

  treerev=(res *)malloc(node*sizeof(res));

  outd=0;           //=0,1,5,21,...
  ind=node; //=tn-1,tn-5,tn-21,t-85?,...

  for(l=1,p=0 ; p<node ; l<<=1)
    {
      ind-=l*l;
      for(i=0; i<l*l; i++)
	{
	  p++;
	  for(k=0; k<40; k++)
	    treerev[outd+i].entry[k] = tree[ind+i].entry[k];
	  treerev[outd+i].side = tree[ind+i].side;
	  treerev[outd+i].avgy = tree[ind+i].avgy;
	}
      outd+=l*l;
    }

  curimg.nnodes=node;

  fwrite(&curimg,sizeof(tfile),1,datfile);
  fwrite(treerev,node*sizeof(res),1,datfile);
  printf("%d\n",node);
  fclose(datfile);
  
  free(tree);
  free(treerev);
  free(arr);
  free(scr);
  free(scr2);
  free(line);
}