#include "includes.h"
 
 Parameter *Encode_Parm;
 BitMap3D  *Encode_Block;
 BitMap3D  *Encode_Image;
 unsigned char Encode_DeltaQuant;
 int Encode_PostTrans;
 
 FILE *fp;
 
/**************************************|****************************************
Routine   : GetClassCount
Input  par: FeatureSpace *S (pointer to feature space)
            int *Index (pointer to index in the feature space)
Output par: unsigned long (count of subclass members)
Function  : Finds the number of nodes in a given subclass of a feature space.
***************************************|***************************************/
 
 unsigned long GetClassCount(FeatureSpace *S,int *Index)
  {
   int i;
   Grid *tmp=S->GridArray;
   GridTerm *T;
   for(i=0;i<S->Dimension;i++) tmp=tmp->Next[Index[i]];
   T=(GridTerm *)tmp;
   return T->Count;
  }
 
/**************************************|****************************************
Routine   : FindNearest
Input  par: ListNode *List (pointer to list of domain blocks)
            PoolNode *Range(range block to compare to)
            int dim  (feature dimension)
            int MaxList (maximum nodes in list)
Output par: ListNode * (pointer to found list)
Function  : Finds the nearest domain blocks for a given range block in the 
            n-dimensional feature space.
***************************************|***************************************/
 
 ListNode *FindNearest(ListNode *List,PoolNode *Range,int dim, int MaxList)
  {
   ListNode *s=List,*NewList=GimmeAListNode();/* this routine has to be very quick since its */
   ListNode *Tmp,*p;                          /* called at least once for every edge range block */
   register unsigned int i;
   register unsigned long d;
   unsigned int Found=1;
   register long t;
   
   NewList->Next=NewList->Pred=NewList;
   NewList->Domain=List->Domain;
   
   d=0;   /* find first distance */
   i=dim;
   while(i--) 
    {
     t=Range->Distance[i]-List->Domain->Distance[i];
     if (t>0) d += t; /* Manhattan distance for speed */
     else d -=t;
    }
   
   NewList->L=d; /* initialize new list */
   List=List->Next;
   
   while(List!=NULL)
    {
     d=0;
     i=dim;
     while(i--)
      {
       t=Range->Distance[i]-List->Domain->Distance[i];
       if (t>0) d += t; /* Manhattan distance for speed */
       else d -=t;
      }
     
     p=NewList;
     
     if(Found<=MaxList)
      {
       Found++;
       
       Tmp=GimmeAListNode();
       Tmp->Domain=List->Domain;
       Tmp->L=d;
       
       while((p->L<d) && (p->Next!=NewList)) p=p->Next;
       
       if(p->L<d) p=p->Next; /* check for last in list */
       
       Tmp->Pred=p->Pred;
       Tmp->Pred->Next=Tmp;
       Tmp->Next=p;
       p->Pred=Tmp;
       if (d<=NewList->L) NewList=Tmp; /* new list start */
      }
     else
     if (p->Pred->L>d)
      {
       while(p->L<d) p=p->Next;
       
       Tmp=NewList->Pred; /* last one in list */
       Tmp->Domain=List->Domain;
       Tmp->L=d;
       
       if(p->Next!=NewList && p!=NewList) /* new element not last in list */
        {
         NewList->Pred=Tmp->Pred;     /* reset last list meber */
         NewList->Pred->Next=NewList;
         
         Tmp->Pred=p->Pred;   /* insert member */
         Tmp->Pred->Next=Tmp;
         Tmp->Next=p;
         p->Pred=Tmp;
        }
       if (d<=NewList->L) NewList=Tmp; /* new list start */          
      }
     List=List->Next;
    }
   
   FreeMeAList(s);
   return NewList;
  }
 
/**************************************|****************************************
Routine   : BuildList
Input  par: FeatureSpace *S (pointer to a feature space)
            int *Index (pointer to index in a feature space)
Output par: ListNode * (pointer created list)
Function  : Takes a subclass of a feature space (pointed out by the index) and
            builds a list of the blocks in the class.
***************************************|***************************************/
 
 ListNode *BuildList(FeatureSpace *S,int *Index)
  {
   int i;
   ListNode *Start,*List,*tmpList;
   Grid *tmp=S->GridArray;
   GridTerm *T;
   
   for(i=0;i<S->Dimension;i++) tmp=tmp->Next[Index[i]];
   T=(GridTerm *)tmp;
   
   if (T->Count==0) return (ListNode *)NULL;
   
   tmpList=T->Class;
   Start=List=GimmeAListNode();
   List->Info=T->Count;
   List->Domain=tmpList->Domain;
   tmpList=tmpList->Next;
   
   while (tmpList!=NULL)
    {
     ListNode *pp=GimmeAListNode();
     
     List->Next=pp;
     pp->Pred=List;
     List=List->Next;
     List->Domain=tmpList->Domain;
     tmpList=tmpList->Next;
    }
   
   List->Next=Start;
   Start->Pred=List;
   
   return Start;
  }
 
/**************************************|****************************************
Routine   : InsertList
Input  par: ListNode *List (first list)
            ListNode *Xtra (second list)
Output par: ListNode * (pointer to meged lists)
Function  : Merges to lists.
***************************************|***************************************/
 
 ListNode *InsertList(ListNode *List,ListNode *Xtra)
  {
   if(Xtra==NULL) return List;
   if(List==NULL) return Xtra;
   
   if (Xtra->Info==0) return List;
   if (List->Info==0) return Xtra;
   
   List->Pred->Next=Xtra;
   List->Pred=Xtra->Pred;
   Xtra->Pred->Next=List;
   Xtra->Pred=List->Pred;
   List->Info +=Xtra->Info;
   
   return List;
  }
 
/**************************************|****************************************
Routine   : XtraSearch
Input  par: FeatureSpace *S (pointer the feature space)
            ListNode *List  (list of blocks found so far)
            int *Index (pointer to rangenode index)
            int shell  (the new shell to be search in (in the feature space))
            int dim (feature space dimension)
            int end (end flag)
Output par: ListNode * (pointer to new list of blocks found)
Function  : A really special routine to do an extra search in an n-dimensional
            feature space. In the case where a normal search do not find enough
            blocks this routine is called. It looks searches in the surrounding 
            of the original block index (shell).
***************************************|***************************************/
 
 ListNode *XtraSearch(FeatureSpace *S,ListNode *List,int *Index,int shell,int dim,int end)
  {
   int *NewIndex=GimmeAIntArray(S->Dimension); /* this routine gave me a lot of trouble ! */
   ListNode *L1,*L2;
   register int i,j;
   int d=dim,flag;
   d++;
   
   for(j= -shell;j<= shell;j++)
    {
     for(i=0;i<S->Dimension;i++) NewIndex[i]=Index[i];
     NewIndex[dim] += j;
     
     flag=TRUE;
     for(i=0;i<S->Dimension;i++) if (NewIndex[i]<0 || NewIndex[i]>=S->Size) flag=FALSE;
     
     if ((j== -shell || j== shell || end) && d==S->Dimension && flag)
      {
       L1=BuildList(S,NewIndex);
       List=InsertList(List,L1);
      }
     
     if (d<=S->Dimension && flag)
      {
       if (j== -shell || j== shell || end) L2=XtraSearch(S,NULL,NewIndex,shell,d,TRUE);
       else L2=XtraSearch(S,NULL,NewIndex,shell,d,FALSE);
       List=InsertList(List,L2);
      }
    }
   
   FreeMeAIntArray(NewIndex);
   
   return List;
  }
 
/**************************************|****************************************
Routine   : SearchGrid
Input  par: FeatureSpace *DomainS (pointer to domain feature space)
            PoolNode *RangeNode   (pointer to range block)
            int MinList,int MaxList (list sizes)
Output par: ListNode * (pointer to list of domain blocks)
Function  : Find the nearest domain blocks for a given range block in an n-
            dimensional feature space.
***************************************|***************************************/
 
 ListNode *SearchGrid(FeatureSpace *DomainS,PoolNode *RangeNode,int MinList,int MaxList)
  {
   ListNode *SearchList;
   int Dim=DomainS->Dimension;
   int i,c=0;
   
   SearchList=BuildList(DomainS,RangeNode->Index);
   
   while (((SearchList==NULL) || (SearchList->Info<MinList)) && (c<10))    
    {
     c++; /* he he : c plus plus */
     SearchList=XtraSearch(DomainS,SearchList,RangeNode->Index,c,0,FALSE);
    }
   if (c>1) vprintf(stderr," %d",c);
   
   if (SearchList==NULL)
    {
     vprintf(stderr,"\n search failed: ");
     for(i=0;i<Dim;i++) vprintf(stderr," %d ",RangeNode->Index[i]);
     exit(10);
     /* for(i=0;i<Dim;i++) I[i]=(DomainS->Size)>>1;
     SearchList=ExtraSearch(DomainS,SearchList,I,1);
     vprintf(stderr,".",SearchList->Info);*/
    }
   
   SearchList->Pred->Next=NULL; /* terminate list */
   if (SearchList->Info>MaxList)
   SearchList=FindNearest(SearchList,RangeNode,DomainS->Dimension,MaxList);
   
   return SearchList;
  }
 
/**************************************|****************************************
Routine   : FindTransformation3D
Input  par: int XPos,int YPos, int ZPos (position of block to be encoded)
            int NSquare (the quadtree level)
            int BlockSize (range blocksize)
            PoolStructure *Pool (pool structure !)
            int MinList,int MaxList (list sizes)
Output par: Transformation * (pointer to the found transformation)
Function  : Encodes a specific range block.
***************************************|***************************************/
 
 Transformation *FindTransformation3D(int XPos,int YPos, int ZPos,int NSquare,
                                      int BlockSize,PoolStructure *Pool,
                                      int MinList,int MaxList)
  {
   register unsigned int x,y,z,D=BlockSize<<1;
   Transformation *Trans=GimmeATransformation();
   PoolNode *RNode;
   ListNode *SearchList,*List;
   float *Features;
   
   Trans->BlockSize=BlockSize;
   Trans->D=D;
   
   RNode=Pool->RangeNodes[XPos/BlockSize][YPos/BlockSize][ZPos/BlockSize];
   Features=RNode->Features;
   
   Encode_Block->XSize=Encode_Block->YSize=Encode_Block->ZSize=BlockSize;
   
   if(RNode->Distance==NULL)
    {
     Trans->Type=SHADEBLOCK;  /* find shade transformation */
     Trans->g0  =Features[MEAN];
    }
   else /* find edge transformtion */
    {
     float Alpha,temp;
     int Deltag;
     register long Bestdm=LONG_MAX,dm,Newdm;
     PoolNode *DNode;
     BitMap3D *tempimg;
     int tn,Count=0;
     register unsigned int i,j,k,x,y,z;
     
     Trans->Type=EDGEBLOCK;
     if(Encode_Parm->TPostProcess>=0)
     List=SearchList=SearchGrid(Pool->DomainSpace,RNode,MinList+POSTLISTEXTRA,MaxList+POSTLISTEXTRA);
     else      
     List=SearchList=SearchGrid(Pool->DomainSpace,RNode,MinList,MaxList);
     
     while ((List->Next!=NULL) && (Count<MaxList) && (List->Next!=SearchList))
      {
       DNode=List->Domain;
       if (DNode->Features[VAR]==0.0) ErrorHandler(OUT_OF_RANGE,"Null variance"); /* should not occur ! */

       tempimg=Pool->SampledDomainBitmap3D;
       x=DNode->x>>1;
       y=DNode->y>>1;
       z=DNode->z>>1; 
       
       temp=Features[STDDEV]/DNode->Features[STDDEV]; 
      /* temp=0.0;
       for(i=0;i<Encode_Block->XSize;i++)    
       for(j=0;j<Encode_Block->YSize;j++)
       for(k=0;k<Encode_Block->ZSize;k++) 
       temp += Encode_Block->Map[i][j][k]*tempimg->Map[x+i][y+j][z+k];

       temp/= (i*j*k);
       temp = (temp-Features[MEAN]*DNode->Features[MEAN])/DNode->Features[VAR];*/

       Alpha=MakeFloat(MakeRange(temp,Encode_Parm),Encode_Parm); /* quant alpha */
       Deltag=(int)(Features[MEAN]-DNode->Features[MEAN]*Alpha);
       
       if (Deltag> Encode_DeltaQuant) Deltag= Encode_DeltaQuant;
       if (Deltag<-Encode_DeltaQuant) Deltag=-Encode_DeltaQuant;
       
       for(i=0;i<Encode_Block->XSize;i++)    
       for(j=0;j<Encode_Block->YSize;j++)
       for(k=0;k<Encode_Block->ZSize;k++) 
       Encode_Block->Map[i][j][k]=(tempimg->Map[x+i][y+j][z+k]*Alpha)+Deltag;
       
       dm=LONG_MAX;
       for(i=1;i<=8;i++)
        {
         Newdm=IsoAnddl23D(Encode_Image,XPos,YPos,ZPos,Encode_Block,BlockSize,i);
         if(Newdm<dm)
          {
           dm=Newdm;
           tn=i;
          }
        }
       
       if(dm<Bestdm)
        {
         Bestdm=dm;
         Trans->Alpha=Alpha;
         Trans->Deltag=Deltag;
         Trans->tn=tn-1;
         Trans->Domain=DNode;
         Trans->Distortion=dm;
        }
       Count++;
       List=List->Next;
       
       /* now check the distortion and postprocess if necessary */
       
       if (Count==MaxList && Encode_Parm->TPostProcess>=0)
       if(Encode_Parm->TPostProcess<Trans->Distortion/(BlockSize*BlockSize*BlockSize))      
        {
         MaxList += POSTLISTEXTRA;
         Encode_PostTrans++;
        }
      }
     
     Trans->Distortion /= (BlockSize*BlockSize*BlockSize); /* normalize dl2 */
     FreeMeAList(SearchList);
    }
   
   if(NSquare>0) /* check if there are more sublevels */
    {
     if (Trans->Type==SHADEBLOCK)
      {
       unsigned int B=BlockSize;
       register unsigned int val=Trans->g0;
       register int tmp;
       register unsigned long Distortion=0;
       
       for(x=0;x<B;x++) /* insert distortion */
       for(y=0;y<B;y++) 
       for(z=0;z<B;z++)
        {
         tmp=Encode_Image->Map[XPos+x][YPos+y][ZPos+z]-val;
         Distortion +=tmp*tmp;
        }
       Trans->Distortion=Distortion/(B*B*B);
      }
     
     if(Trans->Distortion>(unsigned long)Encode_Parm->TMainSub) /* check distortion */
      {
       Trans->Type=NOBLOCK; /* unmark main block */
       
       Trans->Sub=GimmeATransArray3D(2,2,2);
       for(x=0;x<2;x++)
       for(y=0;y<2;y++)
       for(z=0;z<2;z++)
       Trans->Sub[x][y][z]= /* encode ALL sub blocks */
       FindTransformation3D(XPos+x*(BlockSize>>1),YPos+y*(BlockSize>>1),
                            ZPos+z*(BlockSize>>1),NSquare-1,BlockSize>>1,
                            Pool->Next,MinList,MaxList);
      }
    }
   
   return Trans;
  }
 
/**************************************|****************************************
Routine   : Encode3D
Input  par: BitMap3D *Image (pointer to the image to encode)
            int NSquare     (number of octree square levels)
            int StartBlockSize (maximum range block size)
            PoolStructure *Pool(pointer to poolstructure)
            Parameter *Pa   (pointer to parameters)
            int MinList,int MaxList (list sizes)
Output par: Transformation **** (pointer to 3d array of transformations)
Function  : Encodes the image. Finds the a fractal transformations for every
            range block.
***************************************|***************************************/
 
 Transformation ****Encode3D(BitMap3D *Image,int NSquare,int StartBlockSize,
                             PoolStructure *Pool,Parameter *Pa,int MinList,int MaxList)
  {
   register unsigned int x,y,z;
   Transformation ****FCCodes=GimmeATransArray3D(Image->XSize/StartBlockSize,
                                                 Image->YSize/StartBlockSize,
                                                 Image->ZSize/StartBlockSize);
   vprintf(stderr,"\nEncoding...");
   
   Encode_Parm=Pa; /* make parmeters global for module */  
   Encode_DeltaQuant=1<<(Pa->DBits-1)-1;
   Encode_Block=GimmeABitMap3D(StartBlockSize,StartBlockSize,StartBlockSize,Image->ImgType);
   Encode_Image=Image;
   Encode_PostTrans=0;
   
   for(x=0;x<Image->XSize/StartBlockSize;x++)
    {
     for(y=0;y<Image->YSize/StartBlockSize;y++)
     for(z=0;z<Image->ZSize/StartBlockSize;z++)
     FCCodes[x][y][z]=FindTransformation3D(x*StartBlockSize,y*StartBlockSize,z*StartBlockSize,
                                           NSquare,StartBlockSize,Pool,MinList,MaxList);
     
     if (Quiet) fprintf(stderr,"%3.1f%c ",(float)(x*100)/Image->XSize,'%');
     else fprintf(stderr,"\n   x:%3d y:%d z:%d =%3.1f%c",x*StartBlockSize,y*StartBlockSize,
                  z*StartBlockSize,(float)(x*StartBlockSize*100)/Image->XSize,'%');
    }
   
   vprintf(stderr,"\n   Post processed %d blocks",Encode_PostTrans);
   Encode_Block->XSize=StartBlockSize;
   Encode_Block->YSize=StartBlockSize;
   Encode_Block->ZSize=StartBlockSize;
   
   FreeMeABitMap3D(Encode_Block);
   
   return FCCodes;
  }