rtkey.cpp

此文件包含了在linux下实现tpr-tree索引的源代码

//   two rectangles from CT to CT + h 
//
double RTkey::intersectArea_n1 (const RTkey& k, RTperiod h) const
{
   RTkey       xtion;
   double      xarea;
   int         i, j, l;
   struct 
   {
      RTtimeStamp  t;
      int          d;    // dimension 
      int          j;    // 0 or 1 - upper or lower bound
      const RTkey* to;   
   }           c, chp[TPR_MAX_DIMS*2];    // changing points
   int         nchp = 0;                  // number of changing points        
   RTtimeStamp tmin = max(CT, max(t1, k.t1));     
   RTtimeStamp tmax = min(CT + h, min(t2, k.t2));
  
   if (Settings.GetTreeType() == RTsettings::tRtree)
   {   
      float beg, end, kbeg, kend;

      if (tmin > tmax) return 0; 
      xarea = tmax - tmin;   
      for (int i = 0; i < Settings.GetDims(); i++)
      {
         if (Speed(i, 0) == 0.0)
	      {
            beg = coords[i][0];
            end = coords[i][1];
         }
         else if (Speed(i, 0) > 0.0)
              {
                 beg = coords[i][0];
                 end = Coordtp(i, 1, t2);
              }
              else
              {
                 beg = Coordtp(i, 0, t2);
                 end = coords[i][1];
              }
         if (k.Speed(i, 0) == 0.0)
         {
            kbeg = k.coords[i][0];
            kend = k.coords[i][1];
         }
         else if (k.Speed(i, 0) > 0.0)
              {
                 kbeg = k.coords[i][0];
                 kend = k.Coordtp(i, 1, t2);
              }
              else
              {
                 kbeg = k.Coordtp(i, 0, t2);
                 kend = k.coords[i][1];
              }
         if (beg > kend) return 0; 
         if (end < kbeg) return 0;
         xarea *= min(end, kend) - max(beg, kbeg);   
      }   
      return xarea;
   }
   
   assert(!t1 && !k.t1);
   if (!overlapInt_n1a(k, tmin, tmax)) return 0;
      
   assert(t1 == k.t1);
   xtion.t1 = t1;

   // find bounds of xtion in all dimensions and their changing points
   for (i = 0; i < Settings.GetDims(); i++)     
      for (j = 0; j < 2; j++)
      {
	 if ( j && Coordtp(i,j,tmin) < k.Coordtp(i,j,tmin) ||
             !j && Coordtp(i,j,tmin) > k.Coordtp(i,j,tmin))
         {
            xtion.coords[i][j] = coords[i][j]; 
            xtion.Speed(i,j)   = Speed(i,j); 
            if ( j && Coordtp(i,j,tmax) <= k.Coordtp(i,j,tmax) ||
                !j && Coordtp(i,j,tmax) >= k.Coordtp(i,j,tmax)) continue;
            c.to = &k;  
         } 
         else 
         {
            xtion.coords[i][j] = k.coords[i][j]; 
            xtion.Speed(i,j)   = k.Speed(i,j); 
            if ( j && k.Coordtp(i,j,tmax) <= Coordtp(i,j,tmax) ||
                !j && k.Coordtp(i,j,tmax) >= Coordtp(i,j,tmax)) continue;
            c.to = this;                 
         }
         c.d = i;
         c.j = j;
         c.t = (Coordtp(i, j, 0) - k.Coordtp(i, j, 0)) /
               (k.Speed(i, j) - Speed(i, j)); 

         // keep the list of changing points ordered on time
         for (l = nchp; l && chp[l - 1].t > c.t; l--);
         memmove (chp + l + 1, chp + l, (nchp - l) * sizeof(c));
         chp[nchp++] = c;
      }
   
   // Add up the total integral from pieces
   for (xarea = 0, l = 0; l < nchp; l++)
   {
      if (chp[l].t > tmin) xarea += xtion.areaInt_n1(tmin, chp[l].t);
      xtion.coords[chp[l].d][chp[l].j] = chp[l].to->coords[chp[l].d][chp[l].j];
      xtion.Speed(chp[l].d, chp[l].j)  = chp[l].to->Speed(chp[l].d, chp[l].j);
      tmin = chp[l].t;
   }
   // ... and the last piece 
   if (tmax > tmin) xarea += xtion.areaInt_n1(tmin, tmax);
   
   return xarea;    
}  

//-----------------------------------------------------------
RTkey* RTkey::expand (const RTkey& k, RTperiod h) const
{
   RTkey *u  = (RTkey*) k.Copy();
      
   u->t2 = TPR_TS_INF;
   if (Settings.GetTreeType() == RTsettings::tRtree) 
   {
      float beg, end, kbeg, kend;
      for (int i = 0; i < Settings.GetDims(); i++)
      {
         if (Speed(i, 0) == 0.0)
         {
            beg = coords[i][0];
            end = coords[i][1];
         }
         else if (Speed(i, 0) > 0.0)
         {
            beg = coords[i][0];            
            end = Coordtp_br(i, 1, t2);
         }
         else
         {
            beg = Coordtp_br(i, 0, t2);
            end = coords[i][1];
         }
         if (k.Speed(i, 0) == 0.0)
         {
            kbeg = k.coords[i][0];
            kend = k.coords[i][1];
         }
         else if (k.Speed(i, 0) > 0.0)
         {
            kbeg = k.coords[i][0];
            kend = k.Coordtp_br(i, 1, t2);
         }
         else
         {
            kbeg = k.Coordtp_br(i, 0, t2);
            kend = k.coords[i][1];
         }
         
         u->coords[i][0] = min(beg, kbeg); 
         u->coords[i][1] = max(end, kend);
         u->Speed(i, 0) = u->Speed(i, 1) = 0; 
      }
      u->t1 = min(t1, k.t1); 
      u->t2 = max(t2, k.t2); 
   }
   else 
      if (Settings.IsCTUnionMode())
      {
         u->t1 = CT; 
         for (int i = 0; i < Settings.GetDims(); i++)
         {
            u->coords[i][0] = min(k.Coordtp_br(i, 0, CT), Coordtp_br(i, 0, CT));
            u->coords[i][1] = max(k.Coordtp_br(i, 1, CT), Coordtp_br(i, 1, CT));
            u->Speed(i, 0)  = min(k.Speed(i, 0),  Speed(i, 0));
            u->Speed(i, 1)  = max(k.Speed(i, 1),  Speed(i, 1));
         }   
         u->SetRefTS_br(t1);
      }
      else
         for (int i = 0; i < Settings.GetDims() * 2; i++)
         {
            u->coords[i][0] = min(k.coords[i][0], coords[i][0]);
            u->coords[i][1] = max(k.coords[i][1], coords[i][1]);
         }   
   
   return u;
}

//-----------------------------------------------------------
double RTkey::area_n (RTperiod d) const
{
   double area = 1;

   if (expired()) return 0;   
   for (int i = 0; i < Settings.GetDims(); i++)
      area *= CoordT(i, 1, d) - CoordT(i, 0, d);

   return area;
}

//-----------------------------------------------------------
//   Returns the integral of area from CT to CT + h
// 
double RTkey::area_n1  (RTperiod h) const
{
   RTtimeStamp tmin, tmax;
   
   if (Settings.GetTreeType() == RTsettings::tRtree) { tmin = t1; tmax = t2; }
   else                  
   { 
      tmin = CT; 
      tmax = min(CT + h, t2); 
   }
    
   return tmin < tmax ? areaInt_n1 (tmin, tmax) : 0; 
}

//-----------------------------------------------------------
double RTkey::areaInt_n1 (RTtimeStamp tmin, RTtimeStamp tmax) const
{
   double dx1, dx2, dx3;
   double dv1, dv2, dv3;
   double H = tmax - tmin;

   if (Settings.GetTreeType() == RTsettings::tRtree)
   {
      dx2 = dx3 = 1;
      switch(Settings.GetDims())
      {
      case 3: 
         dx3 = coords[2][1] - coords[2][0]; 
      case 2: 
         dx2 = coords[1][1] - coords[1][0]; 
      case 1: 
         dx1 = coords[0][1] - coords[0][0]; 
      }   
      return H*dx1*dx2*dx3;
   }
   switch(Settings.GetDims())
   {
   case 3: 
      dx3 = Coordtp(2, 1, tmin) - Coordtp(2, 0, tmin); 
      dv3 = Speed (2, 1) - Speed (2, 0); 
   case 2: 
      dx2 = Coordtp(1, 1, tmin) - Coordtp(1, 0, tmin); 
      dv2 = Speed (1, 1) - Speed (1, 0); 
   case 1: 
      dx1 = Coordtp(0, 1, tmin) - Coordtp(0, 0, tmin);
      dv1 = Speed (0, 1) - Speed (0, 0); 
   }   
   switch(Settings.GetDims())
   {
   case 1: 
      return H*dx1 + H*H*dv1 / 2;
   case 2: 
      return H*dx1*dx2 + H*H * (dx1*dv2 + dv1*dx2) / 2 + H*H*H*dv1*dv2 / 3;
   case 3: 
      return H*dx1*dx2*dx3 + H*H * (dx1*dx2*dv3 + (dx1*dv2 + dv1*dx2)*dx3) / 2 + 
             H*H*H * ((dx1*dv2 + dv1*dx2)*dv3 + dv1*dv2*dx3) / 3 + H*H*H*H*dv1*dv2*dv3 / 4;
   default: 
      assert(0);
   }      
}

//-----------------------------------------------------------
double RTkey::margin_n(RTperiod d) const
{
   double dx1, dx2, dx3;
 
   if (expired()) return 0;   
   switch(Settings.GetDims())
   {
   case 3: 
      dx3 = CoordT(2, 1, d) - CoordT(2, 0, d); 
   case 2: 
      dx2 = CoordT(1, 1, d) - CoordT(1, 0, d); 
   case 1: 
      dx1 = CoordT(0, 1, d) - CoordT(0, 0, d);
   }   
   switch(Settings.GetDims())
   {
   case 1: 
      return 0;   // has no sence in 1D case 
   case 2: 
      return dx1 + dx2;
   case 3: 
      return dx1 + dx2 + dx3 + dx1*dx2 + dx1*dx3 + dx2*dx3;
   default: 
      assert(0);
   }      
}

//-----------------------------------------------------------
//   Returns the integral of margin from CT to CT + h
// 
double RTkey::margin_n1(RTperiod h) const
{
   double dx1, dx2, dx3;
   double dv1, dv2, dv3;
   double H = (double) min(t2 - CT, h);
   
   switch(Settings.GetDims())
   {
   case 3: 
      dx3 = CoordT(2, 1) - CoordT(2, 0); 
      dv3 = Speed (2, 1) - Speed (2, 0); 
   case 2: 
      dx2 = CoordT(1, 1) - CoordT(1, 0); 
      dv2 = Speed (1, 1) - Speed (1, 0); 
   case 1: 
      dx1 = CoordT(0, 1) - CoordT(0, 0);
      dv1 = Speed (0, 1) - Speed (0, 0); 
   }   
   if (Settings.GetTreeType() == RTsettings::tRtree) 
   {
      H = t2 - t1;
      switch(Settings.GetDims())
      {
      case 1: 
         return dx1 + H ; 
      case 2: 
         return dx1 + dx2 + H + dx1*dx2 + dx1*H + dx2*H; 
      case 3: 
         return dx1 + dx2 + dx3 + H + 
               dx1*dx2 + dx1*dx3 + dx1*H + dx2*dx3 + dx2*H + dx3*H + 
               dx1*dx2*dx3 + dx1*dx2*H + dx1*dx3*H + dx2*dx3*H;
      default: 
        assert(0);
      }
   }
   else   
      switch(Settings.GetDims())
      {
      case 1: 
         return 0;   // has no sence in 1D case 
      case 2: 
         return H * (dx1 + dx2) + H*H * (dv1 + dv2) / 2;
      case 3: 
         return H * (dx1 + dx2 + dx3 + dx1*dx2 + dx1*dx3 + dx2*dx3) + 
                H*H * (dv1 + dv2 + dv3 + dx1*dv2 + dv1*dx2 + dx1*dv3 + 
                       dv1*dx3 + dx2*dv3 + dv2*dx3) / 2 + 
                H*H*H * (dv1*dv2 + dv1*dv3 + dv2*dv3) / 3;
      default: 
        assert(0);
      }      
}

//-----------------------------------------------------------
//   RTkey::dist may be used only in decision making process
//   where distances are only compared to each other, so there is 
//   no need to take square root 
//
double RTkey::dist_n (const RTkey& k, RTperiod d) const
{
   double      sum = 0; 

   for (int i = 0; i < Settings.GetDims(); i++)
      sum += pow((k.CoordT(i, 0, d) + k.CoordT(i, 1, d)) / 2 - 
                 (CoordT(i, 0, d) + CoordT(i, 1, d)) / 2, 2);
   return sum;   
}

//-----------------------------------------------------------
//   RTkey::dist_n1 returns the integral of distance between the 
//   centers of objects from CT to CT + h. True euclidean distance 
//   (with square root) is used, because this changes the ranking of objects,
//   when integrated distance is used. 
//
double RTkey::dist_n1 (const RTkey& k, RTperiod h) const
{
   double dx[TPR_MAX_DIMS], dv[TPR_MAX_DIMS];
   double a, b, c, f, l, m, n;
   int    i;

   if (Settings.GetTreeType() == RTsettings::tRtree)
   { 
      double      sum = pow((k.t1 + k.t2) / 2 - (t1 + t2) / 2, 2); 
      float       beg, end, kbeg, kend;

      for (int i = 0; i < Settings.GetDims(); i++)
      {
         if (Speed(i, 0) == 0.0)
   	   {
            beg = coords[i][0];
            end = coords[i][1];
         }
         else if (Speed(i, 0) > 0.0)
              {
                 beg = coords[i][0];
                 end = Coordtp(i, 1, t2);
              }
              else
              {
                 beg = Coordtp(i, 0, t2);
                 end = coords[i][1];
              }
         if (k.Speed(i, 0) == 0.0)
         {
            kbeg = k.coords[i][0];
            kend = k.coords[i][1];
         }
         else if (k.Speed(i, 0) > 0.0)
              {
                 kbeg = k.coords[i][0];
                 kend = k.Coordtp(i, 1, t2);
              }
              else
              {
                 kbeg = k.Coordtp(i, 0, t2);
                 kend = k.coords[i][1];
              }    
         sum += pow((kbeg + kend) / 2 - (beg + end) / 2, 2);
      }
      return sum;   
   } 

   h = min(min(t2 - CT, k.t2 - CT), h);
   for (i = 0; i < Settings.GetDims(); i++)
      dx[i] = (k.CoordT(i, 0) + k.CoordT(i, 1)) / 2 - (CoordT(i, 0) + CoordT(i, 1)) / 2;
   for (i = 0; i < Settings.GetDims(); i++)
      dv[i] = (k.Speed(i, 0) + k.Speed(i, 1)) / 2 - (Speed(i, 0) + Speed(i, 1)) / 2;
  
   for (a = 0, i = 0; i < Settings.GetDims(); i++) a += dv[i] * dv[i];
   for (b = 0, i = 0; i < Settings.GetDims(); i++) b += 2 * dx[i] * dv[i]; 
   for (c = 0, i = 0; i < Settings.GetDims(); i++) c += dx[i] * dx[i];

   if (a == 0.0 && c == 0.0) return 0; 
   if (a == 0.0)             return h * sqrt(c);
   if (c == 0.0)             return h*h * sqrt(a) / 2;  
  
   f = sqrt(a*h*h + b*h + c);
   l = 2*a*h + b;
   m = 4*a*c - b*b;
   n = 2*sqrt(a);

   return (l*f + log(l / n + f) * m / n - b*sqrt(c) - log(b / n + sqrt(c)) * m / n) / (4*a);
}

#ifdef PRINTING_OBJECTS
void RTkey::Print(ostream& os) const
{
   int i;
   bool point = true;

   // awkward way to check if we have a point and not a rectangle
   for (int i = 0; i < Settings.GetDims()*2 && point; i++)
      if (coords[i][0] != coords[i][1]) point = false;

   os << '"' << '(' << t1 << ';';
   if (point) 
   { 
      for (i = 0; i < Settings.GetDims(); i++) 
         os << (i ? "," : "") << coords[i][0];
      os << ';';
      for (i = 0; i < Settings.GetDims(); i++) 
	      os << (i ? "," : "") << Speed(i, 0);
   }
   else
   {
      for (i = 0; i < Settings.GetDims(); i++) 
         os << (i ? "," : "") << coords[i][0] << '~' << coords[i][1];
      os << ';';
      for (i = 0; i < Settings.GetDims(); i++) 
         os << (i ? "," : "") << Speed(i, 0) 
            << '~' << Speed(i, 1);
   }  
   os << ")\"";
}
#endif