Changeset 14238

Timestamp:

01/11/13 14:16:33 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: testing faster closest point search

Location:

issm/trunk-jpl/src/c/Container

Files:

• DataSet.cpp (modified) (1 diff)
• Observations.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/Container/DataSet.cpp

@@ -178 +178 @@
         /*Check index in debugging mode*/
         _assert_(this!=NULL);
+        _assert_(offset>=0);
         _assert_(offset<this->Size());
 

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -129 +129 @@
         Observation *observation  = NULL;
 
-        /*If radius is not provided or is 0, return all observations*/
+        /*First, find closest point in Quadtree (fast but might not be the true closest obs)*/
+        this->quadtree->ClosestObs(&index,x_interp,y_interp);
+        _assert_(index>=0);
+        observation=dynamic_cast<Observation*>(this->GetObjectByOffset(index));
+        hmin2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
+
+        /*Compute radius square*/
         if(radius==0) radius=this->quadtree->root->length;
-
-        /*Compute radius square*/
         radius2 = radius*radius;
+        if(hmin2<radius2){
+                radius =sqrt(hmin2);
+                radius2=hmin2;
+        }
 
         /*Find all observations that are in radius*/
@@ -140 +148 @@
                 observation=dynamic_cast<Observation*>(this->GetObjectByOffset(indices[i]));
                 h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
-
-                if(i==0){
+                if(h2<hmin2){
                         hmin2 = h2;
                         index = i;
-                }
-                else{
-                        if(h2<hmin2){
-                                hmin2 = h2;
-                                index = i;
-                        }
                 }
         }  
@@ -175 +176 @@
         for(int i=0;i<n;i++){
                 this->ClosestObservation(&xi,&yi,&obs,x[i],y[i],radius);
-                distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
+                if(xi==UNDEF && yi==UNDEF)
+                 distances[i]=UNDEF;
+                else
+                 distances[i]=sqrt( (x[i]-xi)*(x[i]-xi) + (y[i]-yi)*(y[i]-yi) );
         }
 }/*}}}*/