drawocclusion.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- MegaWave2 Command -----------------------------*/
/* mwcommand
name = {drawocclusion};
version = {"1.21"};
author = {"Simon Masnou"};
function = {"Interactive creation of occlusions on an image"};
usage = {
   'g':[gray=255]->gray         "Occlusions gray level in the occluded image (default: 255)",
   'h':hole_image<-hole_image   "Occluded image",
   'z':zoom->zoom     "zoom factor for image display (even integer)",
   cimage->image      "Input image",
   label<-labelimage  "Output image of labelled occlusions",
   notused->window    "Window for image viewing (internal use)"
 };
*/
/*----------------------------------------------------------------------
 v1.21: revision, fixed mw_delete_point_curve and czoom bugs (S.Masnou)
----------------------------------------------------------------------*/


#include <stdio.h>
#include "mw.h"

/* We need the windows facility */
#include "window.h"

/**********************************************************************/
/* IMPORTANT NOTE : at the end of the labeling process, labels of each 
   connected component are sorted increasingly according to their
   order of appearance */
/**********************************************************************/

/* We assume than there are no more than LABELS_NUMBER labels obtained at first pass
 (before updating) */
#define LABELS_NUMBER 100000

int u_compar_i(u,v)  /*  Called by function qsort for sorting decreasingly */
int *u,*v;
  {
    if ((*u)<(*v)) return (1);
    if ((*u)==(*v)) return (0);
    return (-1);
  }

/****************************************************************************/
/*
The algorithm used here is the following : let us consider the case where we want
to compute the connected components of the set of points with value g (the principle
is analogous for the computation of the connected components of the complement).
We create a subsidiary image, one pixel larger in each direction than the input 
image, where all pixels of input frame with value g have value 1, while the others 
are set to 0. Then we scan the interior of this subsidiary image (borders are not
taken into account) from top to bottom and left to right (only one scan is enough). 
Each pixel with value 1 whose 4 upper and left pixels (in case of 8-connectivity) 
are 0 is given a new label (as if it was a new component).
If some of these 4 pixels are nonzero then the current pixel is given the greatest label
among these 4 and we update the transcoding table which tells us which label is 
connected to which one, i.e. which components are part of the same connected components.
Once the image has been totally scanned, we simply have to replace the value (label)
of each nonzero pixel by the smallest label to which it is associated.
*/
/****************************************************************************/

void mise_a_jour_transcode(transcode,a,b)

int transcode[LABELS_NUMBER];
int a,b;

  {
    if (!b) return;
    if (transcode[a]==0)
      transcode[a]=b;
    else
      {
	if (b!=transcode[a])
	  if (b>transcode[a])
	    mise_a_jour_transcode(transcode,b,transcode[a]);
	  else
	    {
	      mise_a_jour_transcode(transcode,transcode[a],b);
	      transcode[a]=b;
	    }
      }
  }
    
/****************************************************************************/

/* In the transcoding table, a label l1 may be associated with another label l2,
itself associated with a third label l3, etc... The following function simply
performs a recursive association of each label l1,l2,l3 with the smallest possible
label in the chain. */

void refresh(transcode,refresh_transcode,i,first)

int transcode[LABELS_NUMBER],refresh_transcode[LABELS_NUMBER];
int i;
int *first;  
    
  {
    if (transcode[i])
      {
	refresh(transcode,refresh_transcode,transcode[i],first);
	refresh_transcode[i]=*first;
      }
    else
      *first=i;
  }

/****************************************************************************/

void fconnected(In,line,col,FOREGROUND,NUMBER,not_writing_on,complement,connectivity)
float *In;
int line,col;
float FOREGROUND;
int *NUMBER;
char *not_writing_on,*complement,*connectivity;

  {
    int Line=line+2,Col=col+2;
    int *Output;
    register float *ptrin;
    register int *ptrout;
    register int dx,dy,i;
    int kernel[4];
    int label,first;
    int transcode[LABELS_NUMBER]; /* transcodage table obtained at first pass and constructed 
				     following the rule : if a point is connected with two different labels 
				     (more than two is impossible) one writes in the table the link
				     label_max->label_min if label_max is not yet linked, else, say we
				     already have label_max->label' one makes by recursivity 
				     label'->label_min or label_min->label' depending on 
				     max(label',label_min)
				     */
    int refresh_transcode[LABELS_NUMBER];  /* same as transcode except that each label is linked with
					      the smallest value in the chain to which it belongs */
    int normalize_transcode[LABELS_NUMBER];
    register int *ptr_trans,*ptr_ref,*ptr_norm;
    int norme;

    /* The frame used for processing is 2 pixels wider and higher than the input to avoid
       problems with borders */

    Output=(int*)malloc(Line*Col*sizeof(int));
    if (Output==NULL)
      mwerror(FATAL,1,"Not enough memory !\n");
    
    for (i=0,ptr_trans=transcode,ptr_ref=refresh_transcode,
	 ptr_norm=normalize_transcode;i<LABELS_NUMBER;i++,ptr_trans++,ptr_ref++,ptr_norm++)
      {*ptr_trans=0;*ptr_ref=0;*ptr_norm=0;}
    
    for (dy=0,ptrout=Output;dy<Col;dy++,ptrout++) *ptrout=0; 
    for (dx=2,ptrin=In;dx<Line;dx++)
      {
	*ptrout=0;ptrout++;
	for (dy=2;dy<Col;dy++,ptrout++,ptrin++)
	  if (complement)
	    if (*ptrin!=FOREGROUND) *ptrout=1;
	    else *ptrout=0;
	  else
	    if (*ptrin==FOREGROUND) *ptrout=1;
	    else *ptrout=0;
	*ptrout=0;ptrout++;
      }
    for (dy=0,ptrout=Output;dy<Col;dy++,ptrout++) *ptrout=0; 
   
    label=1;
    for (dx=2,ptrout=Output+(Col+1);dx<Line;dx++,ptrout+=2)
      for (dy=2;dy<Col;dy++,ptrout++)
	if (*ptrout)
	  {
	    if (label>LABELS_NUMBER)
	      mwerror(FATAL,1,"There are more than %d labels !\n",LABELS_NUMBER);
	    /* 'kernel' contains the four upper nearest neighbours of current point
	               * * *
		       * + x
		       x x x
	       if 4-connectivity is used then only 2 neighbors are taken into account */
	    kernel[0]=*(ptrout-1);
	    if (connectivity) kernel[1]=(-1);
	    else kernel[1]=*(ptrout-(Col+1));
	    kernel[2]=*(ptrout-Col);
	    if (connectivity) kernel[3]=(-1);
	    else kernel[3]=*(ptrout-(Col-1));
	    qsort(kernel,4,sizeof(int),u_compar_i);
	    /* kernel values are sorted decreasingly */
	    if (kernel[0]==0)
	      *ptrout=label++;   /* This point is given a new label */
	    else
	      {
		*ptrout=kernel[0];
		if (kernel[1])
		  if (kernel[1]!=kernel[0])
		    if (kernel[1]!=transcode[kernel[0]])
		      mise_a_jour_transcode(transcode,kernel[0],kernel[1]);
	      }
	  }
    first=0;
    for (i=label-1;i>0;i--)
      if ((transcode[i]) && (!(refresh_transcode[i])))
	refresh(transcode,refresh_transcode,i,&first);
    
    
    /* We want each connected component to be labelled between 1 and the number
       of components. This is done in the following step */
    norme=1;
    for (i=1,ptr_ref=refresh_transcode+1;i<label;i++,ptr_ref++){
      if (!(*ptr_ref))
	*ptr_ref=i;
      if (normalize_transcode[*ptr_ref]==0)
	{normalize_transcode[*ptr_ref]=norme++;}}
    *NUMBER=norme-1;

     /* Writing on image */
    if (!not_writing_on)
      {
	for (dx=2,ptrout=Output+(Col+1),ptrin=In;dx<Line;dx++,ptrout+=2)
	  for (dy=2;dy<Col;dy++,ptrout++,ptrin++)
	    {
	      if (*ptrout)
		*ptrin=(float)(normalize_transcode[refresh_transcode[*ptrout]]);
	      else
		*ptrin=0;
	    }
      }
    free(Output);
  }  

/* Number of channels for the polygon: 1 (gray-level) */
#define NB_OF_CHANNELS 1

#define min(A,B)     (((A)>(B)) ? (B) : (A))
#define max(A,B)     (((A)>(B)) ? (A) : (B))
#define Labels_Number 10000

Polygons poly;
Polygon pl,npl;
Point_curve pc,npc;
int max_number_of_points_in_poly=0;

int x0,y0,x1,y1,oldx1,oldy1,dx,dy;
char first_point=0;
char print_mode=0;

/* Return the number of points in a polygon */

int give_number_of_point_in_poly(p)

Polygon p;

{
  Point_curve c;
  int n;

  n=0;
  for (c=p->first; c; c=c->next) n++;
  return(n);
}

void DrawLinePoly(ImageWindow,a0,b0,a1,b1)


Wframe *ImageWindow;
int a0,b0,a1,b1;

{
  WSetTypePencil(W_COPY);

  WSetSpecialColorPencil(ImageWindow);

  if (abs(a0-a1) >= abs(b0-b1))
    {
      WDrawLine(ImageWindow,a0,b0+1,a1,b1+1);	  
      WDrawLine(ImageWindow,a0,b0-1,a1,b1-1);	  
    }
  else
    {
      WDrawLine(ImageWindow,a0+1,b0,a1+1,b1);	  
      WDrawLine(ImageWindow,a0-1,b0,a1-1,b1);	  
    }


  WSetColorPencil(ImageWindow,127); 
  WDrawLine(ImageWindow,a0,b0,a1,b1);  

  WSetColorPencil(ImageWindow,0);
}

void List_of_Poly()

{
  Polygon p;
  Point_curve c;
  int i;

  printf("\n\t\tList of recorded polygons:\n");

  for (p=poly->first,i =1; p; p=p->next, i++)
    {
      printf("\tPoly #%d\n",i);
      for (c=p->first; c; c=c->next)
	printf("\t\t(%d,%d)\n",c->x,c->y);
    }
  printf("\n");
}

void Unzoom_Poly(zoom)
int zoom;
{
  Polygon p;
  Point_curve c;

  for (p=poly->first; p; p=p->next)
    {
      for (c=p->first; c; c=c->next)
	{
	  c->x=c->x/zoom;;
	  c->y=c->y/zoom;
	}
    }
}  


void Help()

{
  printf("\n\t\tHelp on line\n");

  printf("\nMouse:\n");
  printf("\tLeft button: Create a new point in the current coordinates (x,y).\n");
  printf("\tMiddle button: Remove the last recorded point.\n");
  printf("\tRight button: Terminate the polygon. Go ready for the next one.\n");

  printf("\nKeyboard:\n");
  printf("\tQ: Soft Quit (save the polygons).\n");
  printf("\tH: Help.\n");
  printf("\t<SPACE>: List of recorded polygons.\n");
  printf("\tP: Resumes or stops printing point coordinates.\n");
  printf("\tR: Redraw the window according to the recorded polygons.\n");
}


void Redraw(ImageWindow)

Wframe *ImageWindow;

{
  Polygon p;
  Point_curve c;
  int a0,b0,a1,b1;

  WSetTypePencil(W_COPY);
  WRestoreImageWindow(ImageWindow,0,0,dx,dy); 

  for (p=poly->first; p; p=p->next)
      for (c=p->first; c && (c->x != -1); c=c->next)
	{
	  if (c == p->first) {a0 = c->x; b0 = c->y;}
	  else 
	    {