chain2_compute.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'dwtrans2d' 2.1                    */
/*                                                                          */
/*      Copyright (C) 1998-2002  E.Bacry, J.Fraleu, J.Kalifa, E. Le Pennec, */
/*                         W.L. Hwang , S.Mallat, S.Zhong                   */
/*      emails : [e-mail]                             */
/*               [e-mail]                               */
/*               [e-mail]                               */
/*               [e-mail]                             */
/*               [e-mail]                               */
/*               [e-mail]                                   */
/*               [e-mail]                               */
/*                                                                          */
/*..........................................................................*/
/*                                                                          */
/*      This program is a free software, you can redistribute it and/or     */
/*      modify it under the terms of the GNU General Public License as      */
/*      published by the Free Software Foundation; either version 2 of the  */
/*      License, or (at your option) any later version                      */
/*                                                                          */
/*      This program is distributed in the hope that it will be useful,     */
/*      but WITHOUT ANY WARRANTY; without even the implied warranty of      */
/*      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       */
/*      GNU General Public License for more details.                        */
/*                                                                          */
/*      You should have received a copy of the GNU General Public License   */
/*      along with this program (in a file named COPYRIGHT);                */
/*      if not, write to the Free Software Foundation, Inc.,                */
/*      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
/*                                                                          */
/*..........................................................................*/


#include "lastwave.h"

#include "extrema2d.h"  
/*
          7  0  1

          6  *  2

          5  4  3
*/

typedef struct W2_neighbor {
  struct ext2 *neighbor[8];
} *W2_NEIGHBOR;


static W2_NEIGHBOR W2_new_point_neighbor(void)
{
  W2_NEIGHBOR neighbor;
  int i;
  if(!(neighbor = (W2_NEIGHBOR)(Malloc(sizeof(struct W2_neighbor)))))
    Errorf ("Mem. alloc for NEIGHBOR failed\n");
  for (i = 0; i < 8 ; i++)
    neighbor->neighbor[i] = NULL;
  return (neighbor);
}


static W2_delete_point_neighbor(W2_NEIGHBOR neighbor)
{
  int i;
  for (i = 0 ; i < 8; i++) 
    neighbor->neighbor[i] = NULL;

  Free(neighbor);
  neighbor = NULL;
}

/******** fill up the neighbors of points in a extlis ************/

static W2_compute_point_neighbor(EXTLIS2 extlis, W2_NEIGHBOR neighbor,int k)
{
  int ncol = extlis->ncol;
  EXT2 *values = extlis->first;
  neighbor->neighbor[0] = values[k - ncol];
  neighbor->neighbor[4] = values[k + ncol];
  neighbor->neighbor[1] = values[k - ncol + 1];
  neighbor->neighbor[5] = values[k + ncol - 1];
  neighbor->neighbor[2] = values[k + 1];
  neighbor->neighbor[6] = values[k - 1];
  neighbor->neighbor[3] = values[k + ncol + 1];
  neighbor->neighbor[7] = values[k - ncol - 1];

}

  
/* link the points in the same level through previous and next */

static void W2_point_pic_link(WTRANS2 wtrans,int level,int uniform)
{
  EXTLIS2 extlis = wtrans->extrep->array[level];
  int nrow = extlis->nrow,
      ncol = extlis->ncol;
  EXT2 *values = extlis->first;
  EXT2 ext,neighbor,coarse, adjacent, adjacent1, adjacent2, adjacent3;
  W2_NEIGHBOR pneighbor;
  int dir,dir1,dir2;
  int i,j,I;
  LWFLOAT measure, measure_tmp;

  /* obtain from coarser */
  for(i = 1,I = ncol; i < nrow - 1; i++,I += ncol)
    for(j = 1; j < ncol - 1; j++)
      if(ext = values[I+j]) {
        if(ext->coarser && ext->coarser->finer && 
	   (ext->coarser->finer == ext)) {
	  coarse = ext->coarser;
          if(coarse->next) /* if coarser->next and its finer is close */
            if(coarse->next->finer) { 
	      neighbor = coarse->next->finer;
              if(abs(neighbor->x - ext->x) < 2 &&
                 abs(neighbor->y - ext->y) < 2)
                ext->next = neighbor;
	    }
          if(coarse->previous)  /* coarser->previous and its finer is close */
            if(coarse->previous->finer) {
	      neighbor = coarse->previous->finer;
              if(abs(neighbor->x - ext->x) < 2 &&
                 abs(neighbor->y - ext->y) < 2) 
                ext->previous = neighbor;
	    }
	}
      }

  /* by gradiant */
  for(i = 1,I = ncol; i < nrow-1; i++,I += ncol)
    for(j = 1; j < ncol - 1; j++)
      if(ext = values[I+j]) { 
	pneighbor = W2_new_point_neighbor();
	W2_compute_point_neighbor(extlis, pneighbor, I+j);

        W2_direction(ext->arg,&dir,&dir1,&dir2);

	/* among the three closed directions, find next */
	if (ext->next) {
	  adjacent = ext->next;
	  if (uniform)  measure = (LWFLOAT)
	    fabs((double)(ext->next->mag - ext->mag));
	  else measure = ext->next->mag;
	}
	else  adjacent = NULL;

	if (dir == dir1) {
	  dir1 = (dir +1)%8;
	  dir2 = (dir +7)%8;
	}

	adjacent1 = pneighbor->neighbor[dir];
	adjacent2 = pneighbor->neighbor[dir1];
	adjacent3 = pneighbor->neighbor[dir2];

	if (adjacent1) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)
		fabs((double)(adjacent1->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent1;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent1;
	      measure = (LWFLOAT) fabs((double)(adjacent->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent1->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent1;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent1;
	      measure = adjacent->mag;
	    }
	}


	if (adjacent2) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)
		fabs((double)(adjacent2->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent2;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent2;
	      measure = (LWFLOAT)fabs((double)(adjacent2->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent2->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent2;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent2;
	      measure = adjacent->mag;
	    }
	}


	if (adjacent3) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)fabs((double)(adjacent3->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent3;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent3;
	      measure = (LWFLOAT)fabs((double)(adjacent3->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent3->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent3;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent3;
	      measure = adjacent3->mag;
	    }
	}
	if (adjacent) ext->next = adjacent;
	else ext->next = NULL;
	      

	/* among the three opposite directions, find previous */
	if (ext->previous) {
	  adjacent = ext->previous;
	  if (uniform)  measure = (LWFLOAT)
	    fabs((double)(ext->previous->mag - ext->mag));
	  else measure = ext->previous->mag;
	}
	else adjacent = NULL;

	if (dir == dir1) {
	  dir1 = (dir +1)%8;
	  dir2 = (dir +7)%8;
	}

	adjacent1 = pneighbor->neighbor[(dir+4)%8];
	adjacent2 = pneighbor->neighbor[(dir1+4)%8];
	adjacent3 = pneighbor->neighbor[(dir2+4)%8];

	if (adjacent1) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)fabs((double)(adjacent1->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent1;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent1;
	      measure = (LWFLOAT)fabs((double)(adjacent1->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent1->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent1;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent1;
	      measure = adjacent1->mag;
	    }
	}


	if (adjacent2) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)fabs((double)(adjacent2->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent2;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent2;
	      measure = (LWFLOAT)fabs((double)(adjacent2->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent2->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent2;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent2;
	      measure = adjacent2->mag;
	    }
	}


	if (adjacent3) {
	  if (uniform) 
	    if (adjacent) {
	      measure_tmp = (LWFLOAT)fabs((double)(adjacent3->mag - ext->mag));
	      if (measure_tmp < measure) {
		adjacent = adjacent3;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent3;
	      measure = (LWFLOAT)fabs((double)(adjacent3->mag - ext->mag));
	    }
	  else
	    if (adjacent) {
	      measure_tmp = adjacent3->mag;
	      if (measure_tmp > measure) {
		adjacent = adjacent3;
		measure = measure_tmp;
	      }
	    }
	    else {
	      adjacent = adjacent3;
	      measure = adjacent3->mag;
	    }
	}
	if (adjacent) ext->previous = adjacent;
	else ext->previous = NULL;

	if(pneighbor) {
	  W2_delete_point_neighbor(pneighbor);
	  pneighbor = NULL;
	}
	
      }

  /** ext->next->previous == ext **/
  /** ext->previous->next == ext **/

    for(i = 0,I = 0; i < nrow; i++,I += ncol)
      for(j = 0; j < ncol; j++)
	if(ext = values[I+j]) {
	  if(neighbor = ext->next) {
	    if(!neighbor->previous)
	      neighbor->previous = ext;
	    else if(neighbor->previous != ext)
	      ext->next = NULL;
	  }
	  if(neighbor = ext->previous) {
	    if(!neighbor->next)
	      neighbor->next = ext;
	    else if(neighbor->next != ext)
	      ext->previous = NULL;
	  }
	}

}

/** compute chainlis from extlis at level **/

static CHAINLIS2 W2_compute_chain(WTRANS2 wtrans, int level)
{
  EXTLIS2 extlis;
  EXT2 *values;
  EXT2 ext;
  CHAINLIS2 chainlis;
  int nrow, ncol;
  CHAIN2 chain;
  int i, j, I;

  extlis = wtrans->extrep->array[level];
  chainlis = wtrans->chainrep->array[level];
  values = (EXT2 *)extlis->first;
  nrow= chainlis->nrow = extlis->nrow;
  ncol= chainlis->ncol = extlis->ncol;

  chain = NewChain2();
  chainlis->first = chain;

  for(i = 1,I = ncol; i < nrow - 1; i++,I += ncol)
    for(j = 1; j < ncol - 1; j++) {
      if(ext = values[I+j])
        if(ext->chain == NULL) {
	  chain->next = NewChain2();
	  chain->next->previous = chain;
	  chain = chain->next;
          chain->first = ext;
          while(ext->previous && (ext->previous->next == ext)) { 
	    /* walk through previous to beginning */
            if(ext->previous == chain->first ||
               ext->previous->chain) {/*no previous is in a chain,no loop*/
	      /* if (ext->previous->chain) {
		chain->head = ext->previous;
	      } */
	      if (ext->previous == chain->first) ext->previous = NULL;
              break;
	    }
            else
              ext = ext->previous;
          }

	  /* if (!(chain->head)) chain->head = ext->previous; */
	  ext->previous = NULL;
	  ext->chain = chain; 

          ext->chain_index = chainlis->size; 
          chain->first = ext;
          chain->size = 1; 
          chain->mag = ext->mag;

          while(ext->next && (ext->next->previous == ext)) { 
	    /* walk through next to ending */
            if(ext->next == chain->first || /* no loop to the header */
               ext->next->chain) { /* no next in a chain */
	      if (ext->next->chain) chain->last = ext->next;
              ext->next = NULL;
              break;
            }else 
	    ext = ext->next;
            ext->chain = chain;
            ext->chain_index = chainlis->size;  
            chain->size++;
            chain->mag += ext->mag;
          }

	  if (!(chain->last)) chain->last = ext->next;
	  ext->next = NULL;

	  /* a chain has completed */
          chain->mag /= chain->size;  
          chainlis->size++;
	}
    }


  /* chainlis is established */
  chain->next = NULL;
  if (chainlis->first->next)  chainlis->first= chainlis->first->next;
                              
  else { chain->first= NULL;}

  W2_collect_point_in_point_pic(extlis, chainlis);

  return(chainlis);
}


/* compute the chainrep from wtrans from coarse to fine */
void ComputeChainrep2(WTRANS2 wtrans,int uniform)
{
  EXTREP2 extrep;
  CHAINREP2 chainrep;

  CHAINLIS2  chainpic;
  int l;

  extrep = wtrans->extrep;
  chainrep = wtrans->chainrep;
 
 
  for(l = extrep->noct; l >= 1; l--) {
    chainpic= wtrans->chainrep->array[l];
    W2_remove_point_pic_link(extrep->array[l]);
    W2_point_pic_link(wtrans,l, uniform);
    
    wtrans->chainrep->array[l] = W2_compute_chain(wtrans,l); 
   
  }
  chainrep->noct = wtrans->noct;
}


void C_ChainCompute2(char ** argv)
  
{
  WTRANS2 wtrans;
  int uniform = YES;

  argv = ParseArgv(argv,tWTRANS2_,NULL,&wtrans,0);

  if (wtrans ==NULL) wtrans= GetWtrans2Cur();

 
  W2_check_point_repr(wtrans->extrep);

  if(wtrans->chainrep->noct == 0) {

    wtrans->chainrep->noct = wtrans->extrep->noct;
  }
  else ClearChainrep2(wtrans);
  
  ComputeChainrep2(wtrans,uniform);
 
}