vpoint.c

image processing including fourier,wavelet,segmentation etc.

     /*** For each angular precision level ... ***/
     for(i=0; i<n_pl; i++) {
       /* Number of orientations */
       ntheta = (int)rint(pow(2.0,(double)(min_pl+i)));
       if (verbose) fprintf(stderr,
	       "\nBuilding data structure for angular precision = pi/%d\n",
	       ntheta);	       

       /* Create the tiling for this precision level */
       M[i]=N;
       T = Tilings[i] = newTiling(ntheta,&(p[i]),&(p_inf[i]),&(M[i]));
       if (verbose) fprintf(stderr,
	       "  Total number of tiles = %d internal + %d external = %d\n",
	       T[0].nx*T[0].ny,T[1].nx*T[1].ny,
	       T[0].nx*T[0].ny+T[1].nx*T[1].ny);
     }


  /*
     Search for meaningful vanishing regions ...
   */

     VPlist[0] = mw_change_flist(NULL,2,0,4);
     VPlist[1] = mw_change_flist(NULL,2,0,4);
    
     searching_masked = 0;
     NVP[0] = NVP[1] = 0;
     itn = 0;
     
     while (1) { /* main loop */

     if (verbose) fprintf(stderr,
			  "\n\n*** Searching for vanishing regions (iteration # %d)\n\n",
			  itn);
     nvp1 = 0;
     /*** For each angular precision level ... ***/
     for(i=0; i<n_pl; i++) { /* for i */

       T = Tilings[i];

       /*** For each segment, update all tiles it meets ***/
       for (N=0,j=0;j<Segments->size;j++) {
	 if (!searching_masked) N++; /* to search non-masked VPs count all segments */
	 if (SIsValid(Segments,j)>0.0) {
	     TAddSegment(Tilings,i,Segments,j,R,X0,Y0);
	     if (searching_masked) N++; /* to search masked VPs count only remaining segments */
	   }
       }

       /*** Compute NFA, meaningfulness for each tile ***/
       B[i]    = binomial_tail(N,p[i]);
       for (ie=0;ie<=1;ie++)
	 for (ix=0;ix<T[ie].nx;ix++)
	   for (iy=0;iy<T[ie].ny;iy++)
	     if (TIsValid(T,ie,ix,iy)) {
	       TMeaning(T,ie,ix,iy) = -log10(M[i]*B[i][SegListSize(TSegs(T,ie,ix,iy))]);
	     } else {
	       TMeaning(T,ie,ix,iy) = TMDLMeaning(T,ie,ix,iy); /* == -Inf by default */
	     }
       /* ATTN: for infinite tiles use p_inf < p !! */
       if (IdentifyOppositeTiles == 1)
	 /* use 2*p_inf as an upper bound of the probability
	    of a line meeting any of both opposite tiles */
	 Binf[i] = binomial_tail(N,2.0*p_inf[i]);
       else
	 Binf[i] = binomial_tail(N,p_inf[i]);
       ie = EXTERIOR_TILING;
       iy = T[ie].ny-1;
	 for (ix=0;ix<T[ie].nx;ix++)
	     if (TIsValid(T,ie,ix,iy)) {
	       TMeaning(T,ie,ix,iy) = -log10(M[i]*Binf[i][SegListSize(TSegs(T,ie,ix,iy))]);
	     } else {
	       TMeaning(T,ie,ix,iy) =  TMDLMeaning(T,ie,ix,iy); /* == -Inf by default */
	     }

       /*** Vanishing points are local maxima of meaningfulness ***/
       nvp = 0;
       for (ie=0;ie<=1;ie++)
	 for (iy=0;iy<T[ie].ny;iy++)
	   for (ix=0;ix<T[ie].nx;ix++) {
	     meaning = TMeaning(T,ie,ix,iy);

	     nvp += 
	     TVP(T,ie,ix,iy) = ( TIsValid(T,ie,ix,iy) &&
				 (meaning>=threshold) &&
				 (meaning> TMeaning2(T,ie,ix+1,iy+1)) &&
				 (meaning> TMeaning2(T,ie,ix+1,iy  )) &&
				 (meaning> TMeaning2(T,ie,ix+1,iy-1)) &&
				 (meaning>=TMeaning2(T,ie,ix-1,iy+1)) &&
				 (meaning>=TMeaning2(T,ie,ix-1,iy  )) &&
				 (meaning>=TMeaning2(T,ie,ix-1,iy-1)) &&
				 (meaning>=TMeaning2(T,ie,ix  ,iy+1)) &&
				 (meaning> TMeaning2(T,ie,ix  ,iy-1)) );

	     if (IdentifyOppositeTiles == 0) {
             /* Treat the special case of pairs of infinite tiles
		in the following way:
		If they are both meaningful with the same nfa,
		then only keep the one with lowest ix (theta)
		coordinate, as meaningful.
	     */
	     if ((ie==EXTERIOR_TILING) &&
		 (iy==(T[ie].ny-1)) &&
		 TVP(T,ie,ix,iy) &&
		 (meaning==TMeaning2(T,ie,ix  ,iy+1)) &&
		 (ix>TOpposite(T,ix,iy)))
	       TVP(T,ie,ix,iy) = 0;
	     }
	   }
       if (verbose) fprintf(stderr,
       "  Found %d maximal meaningful regions at this precision level\n",
	       nvp);

       nvp1+=nvp;
     }

     /*** Vanishing points must be also
	  maximal meaningful over precision levels ***/
     if (verbose) fprintf(stderr,
	     "\nTesting inclusion relationships between precision levels...\n"
	     );	       

     /* for each meaningful tile t(i,ie,ix,iy) */
     for(i=1; i<n_pl; i++) {
       T = Tilings[i];
       for (ie=0;ie<=1;ie++)
	 for (iy=0;iy<T[ie].ny;iy++)
	   for (ix=0;ix<T[ie].nx;ix++)
	     if ((TMeaning2(T,ie,ix,iy)>=threshold) || TMDLVP(T,ie,ix,iy)) {
	       /* find all coarser intersecting tiles t_k(j,IE[k],IX[k],IY[k]) */
	       for(j=0; j<i; j++) {
		 TT = Tilings[j];
		 n = FindIntersection(T,ie,ix,iy,TT,IE,IX,IY);
		 for (k=0; k<n; k++)
		   if ((TMeaning2(TT,IE[k],IX[k],IY[k])>=threshold) || TMDLVP(TT,IE[k],IX[k],IY[k])) {
		     if (TMeaning2(TT,IE[k],IX[k],IY[k])>TMeaning2(T,ie,ix,iy))
		       TVP(T,ie,ix,iy) = 0; /* t(i,ie,ix,iy) is not maximal */
		     else
		       TVP(TT,IE[k],IX[k],IY[k]) = 0; /* t_k(j,IE[k],IX[k],IY[k]) is not maximal */
		   }
	       }
	     }
     }

     if (all) {
       /* All maximal vanishing regions (no MDL) */
       nvp2 = 0;
       for(i=1; i<n_pl; i++) {
	 T = Tilings[i];
	 for (ie=0;ie<=1;ie++)
	   for (iy=0;iy<T[ie].ny;iy++)
	     for (ix=0;ix<T[ie].nx;ix++)
	       if (TVP(T,ie,ix,iy)) {
		 /* Remove t from the list of valid tiles */
		 TVP(T,ie,ix,iy) = TIsValid(T,ie,ix,iy) = 0;
		 /* Add t to the list of MDL meaningful tiles */
		 TMDLVP(T,ie,ix,iy) = 1;
		 TMDLMeaning(T,ie,ix,iy) = TMeaning(T,ie,ix,iy);
		 nvp2++;
		 tuple[0] = (float) i;
		 tuple[1] = (float) ie;
		 tuple[2] = (float) ix;
		 tuple[3] = (float) iy;
		 FlistAdd(VPlist[searching_masked],tuple);
		 TMDLSegs(T,ie,ix,iy) = newSegListRef(TSegs(T,ie,ix,iy));
	       }
       }
       if (verbose)
	 fprintf(stderr,
		 "  %d (out of %d) regions remaining after inclusion test\n",
		 nvp2,nvp1);
       nvp3 = nvp2;

     } else {
     /*** Minimum Description Length ***/
     if (verbose) fprintf(stderr,
	     "\nMinimum Description Length...\n"
	     );
     nvp3 = 0;
     /* search most meaningful tile */
     max_meaning = threshold-1.0;
     for(i=0; i<n_pl; i++) {
       T = Tilings[i];
       for (ie=0;ie<=1;ie++)
	 for (iy=0;iy<T[ie].ny;iy++)
	   for (ix=0;ix<T[ie].nx;ix++)
	     if (TVP(T,ie,ix,iy)) {
	       CSegs = TSegs(T,ie,ix,iy);
	       nsegs = SegListSize(CSegs);
	       meaning = TMeaning(T,ie,ix,iy);
	       if (meaning>max_meaning) {
		 max_i = i;
		 max_ie= ie;
		 max_ix= ix;
		 max_iy= iy;
		 max_meaning=meaning;
	       }
	     }
     }
     /* If meaningful add it to the list of MDL-meaningful VPs */
     if (max_meaning >= threshold) {
       T = Tilings[max_i];
       /* Remove t from the list of valid tiles */
       TVP(T,max_ie,max_ix,max_iy) = TIsValid(T,max_ie,max_ix,max_iy) = 0;
       /* Add t in the list of MDL meaningful tiles */
       TMDLVP(T,max_ie,max_ix,max_iy) = 1;
       TMDLMeaning(T,max_ie,max_ix,max_iy) = max_meaning;
       tuple[0] = (float) max_i;
       tuple[1] = (float) max_ie;
       tuple[2] = (float) max_ix;
       tuple[3] = (float) max_iy;
       FlistAdd(VPlist[searching_masked],tuple);
       /* For each segment j meeting this tile t ... */
       CSegs = TSegs(T,max_ie,max_ix,max_iy);
       for (it=SegListBegin(CSegs); !SegListEnd(CSegs,it); SegListNext(it)) {
	 /* Remove segment j from the list of Segments
	    to be considered in the next iteration */
	 j = SegListValue(it);
	 SIsValid(Segments,j) = 0;
	 /* Add j to MDLSegs of t */
#ifdef HasSTL
	 SegListInsert(TMDLSegs(T,max_ie,max_ix,max_iy),j);
#else
	 tuple[0] = (float) j;
	 FlistAdd(TMDLSegs(T,max_ie,max_ix,max_iy),tuple);
#endif
	   
       }
       nvp3++;
     } /* end if max_meaning >= threshold */

     if (verbose) fprintf(stderr,
	     "  %d (out of %d) regions remaining after MDL.\n",
	     nvp3,nvp1);

     } /* end: if (all) ... else (MDL) ... */

     if (nvp3==0) /* no more (MDL-maximal) meaningful tiles */
       if (masked && !searching_masked)
	 searching_masked = 1; /* look for masked MDL-maximal meaningful tiles */
       else
	 break; /* stop searching */

     /* delete all tile-segment intersections for next iteration */
     for (j=0;j<Segments->size;j++)
       if (SIsValid(Segments,j)>0.0)
	 FlistFlush(STiles(Segments,j));
     for(i=0; i<n_pl; i++) {
       T = Tilings[i];
       for (ie=0;ie<=1;ie++)
	 for (iy=0;iy<T[ie].ny;iy++)
	   for (ix=0;ix<T[ie].nx;ix++)
	     if (TIsValid(T,ie,ix,iy))
	       FlistFlush(TSegs(T,ie,ix,iy));
       free(B[i]);
       free(Binf[i]);
     }

     NVP[searching_masked] += nvp3;

     itn++;

     } /* main loop */


     /*** Prepare output ***/
     output = mw_change_flist(output,NVP[0]+NVP[1],0,10);
     if (segs)
       segs = mw_change_flists(segs,NVP[0]+NVP[1],0);

     for (searching_masked=0; searching_masked<=((masked)?1:0); searching_masked++) {

       if (verbose) {
	 if (all)
	   msg1 = "maximal";
	 else
	   msg1 = "MDL maximal";
	 
	 if (searching_masked)
	   msg2 = "masked ";
	 else
	   msg2 = "";

	 fprintf(stderr,
		 "\n*** Found %d %s%s meaningful vanishing regions:\n",
		 NVP[searching_masked],msg2,msg1);
       }

       for (k=0; k<VPlist[searching_masked]->size; k++) {

	 i  = (int) ((VPlist[searching_masked])->values[k*4+0]);
	 ie = (int) ((VPlist[searching_masked])->values[k*4+1]);
	 ix = (int) ((VPlist[searching_masked])->values[k*4+2]);
 	 iy = (int) ((VPlist[searching_masked])->values[k*4+3]);

	 T = Tilings[i];

	 if (ie == INTERIOR_TILING) {
	   if (verbose) fprintf(stderr,"\nVP%d:\n",output->size);
	   x1 = x4 = (T[ie].x[ix  ]*R + X0);
	   x2 = x3 = (T[ie].x[ix+1]*R + X0);
	   y1 = y2 = (T[ie].y[iy  ]*R + Y0);
	   y3 = y4 = (T[ie].y[iy+1]*R + Y0);
	   meaning = TMeaning(T,ie,ix,iy);
	   output->values[10*(output->size)+0] = x1;
	   output->values[10*(output->size)+1] = y1;
	   output->values[10*(output->size)+2] = x2;
	   output->values[10*(output->size)+3] = y2;
	   output->values[10*(output->size)+4] = x3;
	   output->values[10*(output->size)+5] = y3;
	   output->values[10*(output->size)+6] = x4;
	   output->values[10*(output->size)+7] = y4;
	   output->values[10*(output->size)+8] = meaning;
	   output->values[10*(output->size)+9] = (i+min_pl);
	   output->size++;
	   if (segs)
	     segs->list[segs->size++] = SegListCopy(TMDLSegs(T,ie,ix,iy));

	   if (verbose) fprintf(stderr,
				"   precision level = pi/%g\n",
				pow(2.0,(i+min_pl)));
	   if (verbose) fprintf(stderr,
				"   meaningfulness  = %g\n",
				meaning);
	   if (verbose) fprintf(stderr,
				"   meaningfulness after MDL = %g\n",
				TMDLMeaning(T,ie,ix,iy));
	   if (verbose) fprintf(stderr,
				"   normalized coordinates: x in [%g,%g], y in [%g,%g]\n",
				T[ie].x[ix  ],T[ie].x[ix+1],
				T[ie].y[iy  ],T[ie].y[iy+1]);
	   if (verbose) fprintf(stderr,
				"   pixel coordinates: x in [%g,%g], y in [%g,%g]\n",
				x1,x2,
				y1,y3);
	 } else { /* ie == EXTERIOR_TILING */
	   if (verbose) fprintf(stderr,"\nVP%d:\n",output->size);
	   t1 = t4 = T[ie].x[ix  ];
	   t2 = t3 = T[ie].x[ix+1];
	   d1 = d2 = T[ie].y[iy  ];
	   d3 = d4 = (iy==(T[ie].ny-1))? 2.0*T[ie].y[iy  ] : T[ie].y[iy+1];
	   polar2cart(t1,d1,R,X0,Y0,&x1,&y1);
	   polar2cart(t2,d2,R,X0,Y0,&x2,&y2);
	   polar2cart(t3,d3,R,X0,Y0,&x3,&y3);
	   polar2cart(t4,d4,R,X0,Y0,&x4,&y4);
	   meaning = TMeaning(T,ie,ix,iy);
	   output->values[10*(output->size)+0] = x1;
	   output->values[10*(output->size)+1] = y1;
	   output->values[10*(output->size)+2] = x2;
	   output->values[10*(output->size)+3] = y2;
	   output->values[10*(output->size)+4] = x3;
	   output->values[10*(output->size)+5] = y3;
	   output->values[10*(output->size)+6] = x4;
	   output->values[10*(output->size)+7] = y4;
	   output->values[10*(output->size)+8] = meaning;
	   output->values[10*(output->size)+9] = (i+min_pl);
	   output->size++;
	   if (segs)
	     segs->list[segs->size++] = SegListCopy(TMDLSegs(T,ie,ix,iy));
	   
	   if (verbose) fprintf(stderr,
				"   precision level = pi/%g\n",
				pow(2.0,(i+min_pl)));
	   if (verbose) fprintf(stderr,
				"   meaningfulness  = %g\n",
				meaning);
	   if (verbose) fprintf(stderr,
				"   meaningfulness after MDL = %g\n",
				TMDLMeaning(T,ie,ix,iy));
	   if (verbose) fprintf(stderr,
				"   normalized coordinates: theta in [%g,%g], rho in [%g,%g]\n",
				t1,t2,
				d1,d3);
	   if (verbose) fprintf(stderr,
				"   pixel coordinates: x in [%g,%g], y in [%g,%g]\n",
				x1,x2,
				y1,y3);
	 }
       }
     }


     /* Cleanup */

     deleteSegments(Segments);
     for(i=0; i<n_pl; i++) {
       deleteTiling(Tilings[i]);
       free(B[i]);
       free(Binf[i]);
     }
     free(Tilings);
     free(B);
     free(Binf);
     free(M);

     /* screen output */
     *maskedVPs = NVP[1];
     return NVP[0];


}