www.pudn.com > jseg.rar > segment.c

#include 
#include 
#include 
#include 
#include "mathutil.h"
#include "imgutil.h"
#include "memutil.h"
#include "ioutil.h"
#include "segment.h"

int MM=0;
unsigned char *rmap9;
char tmpfname[200];

int segment(unsigned char *rmap0,unsigned char *cmap,int N,int nt,int ny,int nx,
    unsigned char *RGB,char *outfname,char *exten,int type,int dim,int NSCALE,
    float displayintensity,int verbose,int tt)
{
  float MINRSIZE[5];
  int scale[5],offset[5],step[5],MAXSCALE,MINSCALE,*count,TR,i,j,k,l;
  int oldTR,extraTR,*reg,*reg2,*convert,datasize,autoscale;
  char fname[200];
  short *rmap;

  printf("start segmentation\n");
  scale[0]=32;  offset[0]=2;  step[0]=1;  /* not used */
  scale[1]=64;  offset[1]=4;  step[1]=1;
  scale[2]=128; offset[2]=8;  step[2]=2;
  scale[3]=256; offset[3]=16; step[3]=4;
  scale[4]=512; offset[4]=32; step[4]=8;
  for (i=4;i>=1;i--)
  {
    if (ny*nx>=sqr(scale[i])) { MAXSCALE=i; break; }
  }

  if (i==0) { printf("minimum image size 64x64\n"); return 0; }
  if (nt==1)
  {
    for (i=0;i<=MAXSCALE;i++) MINRSIZE[i]=2.0*sqr(offset[i]);
  }
  if (nt>1) 
  {
    for (i=0;i<=MAXSCALE;i++) MINRSIZE[i]=sqr(offset[i]);
  }

  if (NSCALE>0) 
  {
    autoscale=0;
    MINSCALE=MAXSCALE-NSCALE+1;
    if (MINSCALE<=0) { MINSCALE=1; NSCALE=MAXSCALE-MINSCALE+1; }
  }
  else if (NSCALE==-1)
  {
    autoscale=1;
    if (MAXSCALE>=2)
    {
      MINSCALE = 2;
      NSCALE=MAXSCALE-MINSCALE+1;
      if (NSCALE>2)
      {
        NSCALE = 2;
        MINSCALE=MAXSCALE-NSCALE+1;
      }
    }
    else { MINSCALE = 1; NSCALE = 1; }
  }

  datasize = nt*ny*nx;
  rmap=(short *)calloc(datasize,sizeof(short));
  for (l=0;l=MINSCALE;i--)
  {
    for (l=0;l1) printf("init %d TR=%d %f\n",i,TR,MINRSIZE[i]);

    reg = (int *) calloc(oldTR+1,sizeof(int));
    reg2 = (int *) calloc(TR+1,sizeof(int));
    for (l=0;l1) tempofilt(rmap0,nt,ny,nx,N,cmap,offset,step);
    count = (int *)calloc(TR+1,sizeof(int));
    for (l=0;l1) printf("redo ");
      oldTR=TR;
      convert = (int *)calloc(oldTR+1,sizeof(int));
      for (j=1;j<=oldTR;j++) convert[j] = j;

      for (j=1;j<=oldTR;j++)
      {
        if (count[j]>tt*sqr(scale[i])/8 && reg2[j]==-1) 
        {
          if (verbose || nt>1) printf("j=%d ",j);
          convert[j]=0;
          TR --; 
          for (k=j+1;k<=oldTR;k++) convert[k]--;
        }
      }
      if (TRoldTR-TR)
        {
          reg2 = (int *) realloc(reg2,(TR+extraTR+1)*sizeof(int));
          reg = (int *)calloc(oldTR+1,sizeof(int));
          for (j=1;j<=oldTR;j++) reg2[convert[j]]=reg2[j];
          for (j=TR+1;j<=TR+extraTR;j++) reg2[j]=0; 
        
          for (l=0;l0) 
            {
              k = TR+rmap[l];
              rmap0[l] = k;
              if (reg[j]==0) reg[j] = k; 
              else if (reg[j]==-1) reg2[k]=-1; 
              else if (reg[j]!=k) { reg2[k]=-1; reg2[reg[j]]=-1; reg[j]=-1; }
            }
            else if (rmap[l]==0) rmap0[l]=0;
            else rmap0[l] = convert[j];
          }
          free(reg);
        }
        TR+=extraTR;
      }
      free (convert);

      if (TR>oldTR)
      {
        if (nt>1) tempofilt(rmap0,nt,ny,nx,N,cmap,offset,step);
        count = (int *)realloc(count,(TR+1)*sizeof(int));
        for (j=1;j<=TR;j++) count[j]=0;
        for (l=0;l1) printf("TR=%d \n",TR);
    } while (oldTR!=TR);
    free(count); free(reg2);
    oldTR=TR;
    if (verbose)
      outputEdge(outfname,exten,RGB,rmap0,ny,nx,i,type,dim,displayintensity);
    printf("%d TR=%d\n",i,TR);
  }

  if (autoscale==1 && MINSCALE>1 && NSCALE<3 && nt==1)
  {
    i = MINSCALE-1;
    TR=segment2(rmap0,rmap,i,cmap,N,nt,ny,nx,tt,oldTR,MINRSIZE,offset,step,verbose);
    if (verbose)
    {
      if (verbose || nt>1) printf("TR=%d \n",TR);
      outputEdge(outfname,exten,RGB,rmap0,ny,nx,i,type,dim,displayintensity);
    }
  }

  free(rmap);
  return TR;
}

int segment2(unsigned char *rmap0,short *rmap,int i,unsigned char *cmap,int N,int nt,
    int ny,int nx,int tt,int oldTR,float *MINRSIZE,int *offset,int *step,int verbose)
{
  float oldJ[256],**mapmatrix,*change,**P;
  int TR,j,k,l,*convert,imgsize,*count,extraTR;
  int debug=0;

  imgsize = ny*nx;
  for (k=0;k<256;k++) oldJ[k]=0;
  mapmatrix = (float **)fmatrix(oldTR+1,N+1);
  gettotalJS(cmap,N,ny,nx,rmap0,oldTR,oldJ,mapmatrix,oldTR);

  count = (int *)calloc(oldTR+1,sizeof(int));
  P = (float **)fmatrix(oldTR+1,N);
  for (l=0;l2*MINRSIZE[i])
        change[j]=1.0;
    }
if (debug) 
{
  if (change[j]>1.0 || oldJ[j]>1.0)
  {
    printf("%d ",j);
    for (k=0;k0.01) printf("%1.0f %1.0f; ",mapmatrix[j][k],P[j][k]); 
    }
    printf("%1.0f %4.2f \n",mapmatrix[j][N],oldJ[j]);
  }  
}
  }
/*
printf("%d %d %d\n",rmap0[60*352+140],rmap0[151*352+112],rmap0[202*352+348]);
*/

  TR=oldTR;
  for (j=1;j<=oldTR;j++) convert[j]=j; 
  for (j=1;j<=oldTR;j++)
  {
    if (change[j]>=1.0 || oldJ[j]>1.0) 
    {
      convert[j]=0;
      for (k=j+1;k<=oldTR;k++) convert[k]--;
      TR--;
    }
  }

  for (l=0;l0) rmap[l]=-1; 
    else rmap[l]=0;
  }
  extraTR=segment1(cmap,N,nt,ny,nx,offset,step,rmap,rmap0,oldTR,i,MINRSIZE[i],0,tt);
  for (l=0;l0) rmap0[l]=rmap[l]+TR; 
    else rmap0[l]=convert[rmap0[l]];
  }
  TR += extraTR;
  if (verbose) printf("init %d TR=%d %f\n",i,TR,MINRSIZE[i]);

  free(count);
  free(convert);
  free(change);
  free_fmatrix(P,oldTR+1);
  free_fmatrix(mapmatrix,oldTR+1);
  return TR;
}

int segment1(unsigned char *cmap,int N,int nt,int ny,int nx,int *offset,int *step, 
    short *rmap,unsigned char *rmap0,int oldTR,int i,float MINRSIZE,int redo,int tt)
{
  float *J,*JT,**threshJ1,**threshJ2;
  int j,TR,**done,alldone,datasize,it,imgsize,l;

  imgsize = ny*nx;
  datasize = nt*ny*nx;
  J = (float *)calloc(datasize,sizeof(float));

  for (it=0;it1)
  {
    JT = (float *)calloc(datasize-imgsize,sizeof(float));
    for (it=0;it1) free(JT);

  removehole(rmap,nt,ny,nx,rmap0);
  alldone=getrmap2(rmap,J,nt,ny,nx,TR,oldTR,rmap0,done);
  if (alldone=MAX(i-2,1);j--)
    {
      for (it=0;it0) rmap9[l]=rmap1[l]+1;
    else if (rmap1[l]==0) rmap9[l]=1;
  }
  outputimgraw(tmpfname,rmap9,ny,nx,1);
  free(rmap9);
}
*/
 
  rmap1 += imgsize; rmap0 += imgsize; J += imgsize;
  st=1;
  while (TR==0 && st<=nt-tt)
  {
    TR=getrmap3(rmap1,J,ny,nx,threshJ1[st],0,MINRSIZE,rmap0,n2bgrow,threshJ2[st],
        oldTR,appear[st]);
    rmap1 += imgsize; rmap0 += imgsize; J += imgsize; JT += imgsize;
    st++;
  }
  if (st>nt-tt+1) return TR;
  tracklen = (int *) calloc(TR+1,sizeof(int));
  for (i=1;i<=TR;i++) tracklen[i]=1;

  threshJ3 = (float **)fmatrix(TN,oldTR+1);
  for (k=0;kTR)
            TR4[k]=track1(&(tracklen1[k]),TR,TR2,imgsize,rmap1-imgsize,rmap2[k],
                JT,rmap0,convert[k],&(TR1[k]),&(newTR[k]));

/*
if (MM==1) 
{
  printf("k=%d TR2=%d TR1=%d newTR=%d TR4=%d\n",k,TR2,TR1[k],newTR[k],TR4[k]);
  for (i=1;i<=TR2;i++) printf("%d ",convert[k][i]);
  printf("\n");
}
*/

        }
        piksrtint(TN,TR4,index);
        if (TR4[0]>=0)
        {
          for (k=0;kTR)
            TR4[k]=track1(&(tracklen1[k]),TR,TR2,imgsize,rmap1-imgsize,rmap2[k],
                JT,rmap0,convert[k],&(TR1[k]),&(newTR[k]));
        }

/*
if (MM==1)
{
  rmap9=(unsigned char *)calloc(imgsize,sizeof(unsigned char));
  sprintf(tmpfname,"rmap.%d.%d.gray",it,j);
  for (l=0;l0) rmap9[l]=200;
    else if (rmap0[l]==j) rmap9[l]=100;
  }
  outputimgraw(tmpfname,rmap9,ny,nx,1);
  free(rmap9);  
}
*/

        if (TR1[index[0]]>0)
        {
          if (newTR[index[0]]0) rmap[l]=convert[index[0]][rmap[l]];
            }
            for (l=0;l0) rmap1[l]=convert[index[0]][rmap1[l]];
            }
          }
          for (l=0;l0) 
              rmap1[l]=convert[index[0]][rmap2[index[0]][l]];
          }
          tracklen=(int *) realloc(tracklen,(TR1[index[0]]+1)*sizeof(int));
          for (i=1;i<=TR1[index[0]];i++) tracklen[i]=tracklen1[index[0]][i];
          TR = TR1[index[0]];
if (MM==1)
{
  printf("index=%d TR=%d Track ",index[0],TR);
  for (i=1;i<=TR;i++) printf("%d ",tracklen[i]);
  printf("\n");
}
        }
        for (k=0;k0) rmap9[l]=rmap1[l]+1;
    else if (rmap1[l]==0) rmap9[l]=1;
  }
  outputimgraw(tmpfname,rmap9,ny,nx,1);
  free(rmap9);
}
*/
    rmap1 += imgsize; rmap0 += imgsize; J += imgsize; JT += imgsize;
  }
  free_fmatrix(threshJ3,TN);
  for (k=0;k0) rmap[l]=convert1[rmap[l]];
    }
  }
  free(convert1);
  free(tracklen);

/*
rmap1 = rmap; rmap0 = rmap00;
check1=(int *) calloc(oldTR+1,sizeof(int));
printf("recheck\n");
for (it=0;it0) check1[rmap0[l]]=1;
  }
  for (i=1;i<=oldTR;i++)
  {
    if (appear[it][i] && check1[i]==0) printf("it=%d %d\n",it,i);
  }
  rmap1 += imgsize; rmap0 += imgsize;
}
free(check1);
*/

  free_imatrix(appear,nt); 
  return TR3;
}

int track1(int **tracklen,int TR,int TR2,int imgsize,short *rmap1,short *rmap2,
    float *JT,unsigned char *rmap0,int *convert,int *TR1,int *newTR)
{
  int **neigh,i,j,k,l,l1,m,TR4,*appear,Tappear,*tracked;

/*
printf("TR=%d TR2=%d \n",TR,TR2);
*/

  appear = (int *)calloc(TR+1,sizeof(int));
  neigh = (int **)imatrix(TR+1,TR2+1);
  for (l=0;l0) 
    {
      appear[rmap1[l]]=1;
      if (rmap2[l]>0 && JT[l]<0.2)
      {
        if (rmap0[l]==rmap0[l1]) neigh[rmap1[l]][rmap2[l]] = 1;
      }
    }
  }

  tracked = (int *)calloc(TR+1,sizeof(int));
  TR4=0; Tappear=0;
  for (i=1;i<=TR;i++)
  {
    Tappear+=appear[i];
    for (j=TR+1;j<=TR2;j++)
    {
      if (neigh[i][j] == 1) { (*tracklen)[i]++; TR4++; tracked[i]=1; break; }
    }
  }
  free(appear);
  TR4 = TR4-Tappear;
  *TR1 = TR2;
  for (i=1;i<=TR;i++)
  {
    for (j=TR+1;j<=TR2;j++)
    {
      if (neigh[i][j] == 1)
      {
        for (k=j+1;k<=TR2;k++)
        {
          if (convert[k]>convert[j]) convert[k]--;
        }
        convert[j] = i;
        for (k=i+1;k<=TR;k++)
        {
          if (neigh[k][j] == 1)
          {
            neigh[i][k] = neigh[k][i] = 1;
            neigh[k][j] = 0;
          }
        }
        (*TR1)--;
      }
    }
  }

/*
printf("TR4=%d\n",TR4);
*/
  *tracklen = (int *) realloc(*tracklen,((*TR1)+1)*sizeof(int));
  for (i=TR+1;i<=(*TR1);i++) (*tracklen)[i]=1;
  *newTR=TR;
  do
  {
    m=0;
    for (i=1;i<(*newTR);i++)
    {
      for (j=i+1;j<=(*newTR);j++)
      {
        if (neigh[i][j]==1)
        {
/*
          if (tracked[j]) TR4--;
*/
          TR4--;
          for (k=1;k<=TR2;k++)
          {
            if (convert[k]==j) convert[k]=i;
            else if (convert[k]>j) convert[k]--;
          }
          (*tracklen)[i] = MAX((*tracklen)[i],(*tracklen)[j]);
          for (k=j+1;k<=(*TR1);k++) (*tracklen)[k-1] = (*tracklen)[k];
          for (k=1;k<=(*newTR);k++)
          {
            if (k!=i && k!=j)
            {
              if (neigh[j][k]==1) { neigh[i][k] = neigh[k][i] = 1; }
            }
          }
          for (k=1;k<=(*newTR);k++)
            for (l=j+1;l<=(*newTR);l++) neigh[k][l-1]=neigh[k][l];
          for (k=j+1;k<=(*newTR);k++)
          {
            for (l=1;l<=(*newTR)-1;l++) neigh[k-1][l]=neigh[k][l];
            tracked[k-1]=tracked[k];
          }
          (*newTR)--; (*TR1)--; 
          m = 1;
        }
      }
    }
  } while (m==1);
  free(tracked);
  free_imatrix(neigh,TR+1);
  return TR4;
}

int merge1(unsigned char *rmap,unsigned char *cmap,int N,int nt,int ny,int nx,int TR,
    float threshcolor)
{
  int it,iy,ix,ir,jr,i,mini,minj,newtr,*npt,npttotal;
  float *currentJ,*mergeJ,mindist,**distnpt,**distcolor,**distJ,**unused,**P;
  int loc,loc1,l,datasize,imgsize,*convert,threshtr;
  unsigned char *rmap1,*rmap2;
  float oldoverallJ,overallJ;
int debug=0;

  if (threshcolor>=0) 
  {
    TR=merge(rmap,cmap,N,nt,ny,nx,TR,threshcolor,0);
    return TR;
  }

  threshcolor=0.5; 
  threshtr = TR;

  imgsize = ny*nx;
  datasize = nt*imgsize;
  distnpt=(float **)fmatrix(TR+1,TR+1);
  distJ=(float **)fmatrix(TR+1,TR+1);
  distcolor=(float **)fmatrix(TR+1,TR+1);
  loc=0;
  for (it=0;it1) printf("%d %d %f \n",ir,jr,distJ[ir][jr]);
      }
    }
  }

  convert = (int *)calloc(TR+1,sizeof(int));
  for (i=1;i<=TR;i++) convert[i]=i;
  newtr=TR;
  while (mindist1) 
  printf("min %d %d %f %f %f\n",mini,minj,currentJ[mini],currentJ[minj],mindist);

    overallJ = overallJ-(npt[mini]*currentJ[mini]+npt[minj]*currentJ[minj])/datasize;
    npttotal = npt[mini] + npt[minj];
    for (i=0;i1) 
  printf("%d %d %f %f %f %f\n",mini,ir,distJ[mini][ir],currentJ[mini],currentJ[ir],mergeJ[1]);
        }
      }
    }

    for (i=1;i<=TR;i++) { if (convert[i]>minj) convert[i]--; }

    for (ir=minj+1;ir<=newtr;ir++) 
    {
      npt[ir-1]=npt[ir];
      currentJ[ir-1]=currentJ[ir];
      for (i=0;ithreshtr)
  {
    for (i=1;i<=TR;i++)
    {
      if (convert[i]==minj) convert[i]=mini;
      else if (convert[i]>minj) convert[i]--;
    }

    npttotal = npt[mini] + npt[minj];
    for (i=0;i0)
      {
        if (appear[it-1][rmap[l]] && appear[it+1][rmap[l]]) 
          rmap[l]=rmap[l-imgsize];
      }
    }

    for (l=0;l0)
      {
        if (rmap[l]!=rmap[l-imgsize] && appear[it-1][rmap[l]]) rmap[l]=0;
        else if (rmap[l]!=rmap[l+imgsize] && appear[it+1][rmap[l]]) rmap[l]=0;
      }
    }

    for (l=0;l0)
        {
          if (iy-1>=0) { if (rmap[l-nx]==0) rmap2[l]=0; }
          if (ix-1>=0) { if (rmap[l-1] ==0) rmap2[l]=0; }
          if (ix+1