Merge pull request #33 from leiyangleon/fix_pixel_shift_error

Fix pixel shift error
nasa-jpl · Jun 4, 2021 · 64a3565 · 64a3565
2 parents a2086d7 + e7b17f7
commit 64a3565
Show file tree

Hide file tree

Showing 5 changed files with 19 additions and 12 deletions.
diff --git a/geo_autoRIFT/autoRIFT/__init__.py b/geo_autoRIFT/autoRIFT/__init__.py
@@ -13,4 +13,4 @@
     # this means ISCE support not available. Don't raise error. Allow standalone use
     pass
 
-__version__ = '1.3.0'
+__version__ = '1.3.1'
diff --git a/geo_autoRIFT/autoRIFT/autoRIFT.py b/geo_autoRIFT/autoRIFT/autoRIFT.py
@@ -1372,11 +1372,11 @@ def colfilt(A, kernelSize, option, chunkSize=4):
             A1 = np.lib.pad(A[:,startInds:endInds],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(int((kernelSize[1]-1)/2),int((kernelSize[1]-1)/2))),mode='constant',constant_values=np.nan)
         else:
             if ii == 0:
-                A1 = np.lib.pad(A[:,startInds:endInds+int((kernelSize[1]-1)/2)],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(int((kernelSize[1]-1)/2),0)),mode='constant',constant_values=np.nan)
+                A1 = np.lib.pad(A[:,startInds:np.min((endInds+int((kernelSize[1]-1)/2),A.shape[1]-1))],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(int((kernelSize[1]-1)/2),np.max((0,endInds+int((kernelSize[1]-1)/2)-A.shape[1]+1)))),mode='constant',constant_values=np.nan)
             elif ii == chunkSize-1:
-                A1 = np.lib.pad(A[:,startInds-int((kernelSize[1]-1)/2):endInds],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(0,int((kernelSize[1]-1)/2))),mode='constant',constant_values=np.nan)
+                A1 = np.lib.pad(A[:,np.max((0,startInds-int((kernelSize[1]-1)/2))):endInds],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(np.max((0,0-startInds+int((kernelSize[1]-1)/2))),int((kernelSize[1]-1)/2))),mode='constant',constant_values=np.nan)
             else:
-                A1 = np.lib.pad(A[:,startInds-int((kernelSize[1]-1)/2):endInds+int((kernelSize[1]-1)/2)],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(0,0)),mode='constant',constant_values=np.nan)
+                A1 = np.lib.pad(A[:,np.max((0,startInds-int((kernelSize[1]-1)/2))):np.min((endInds+int((kernelSize[1]-1)/2),A.shape[1]-1))],((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(np.max((0,0-startInds+int((kernelSize[1]-1)/2))),np.max((0,endInds+int((kernelSize[1]-1)/2)-A.shape[1]+1)))),mode='constant',constant_values=np.nan)
 
         B = viewW(A1, kernelSize).reshape(-1,kernelSize[0]*kernelSize[1]).T[:,::1]
 

diff --git a/geo_autoRIFT/geogrid/GeogridOptical.py b/geo_autoRIFT/geogrid/GeogridOptical.py
@@ -138,8 +138,8 @@ def determineBbox(self, zrange=[-200,4000]):
 #        pdb.set_trace()
 
 
-        samples = self.startingX + np.array([0, self.numberOfSamples]) * self.XSize
-        lines = self.startingY + np.array([0, self.numberOfLines]) * self.YSize
+        samples = self.startingX + np.array([0, self.numberOfSamples-1]) * self.XSize
+        lines = self.startingY + np.array([0, self.numberOfLines-1]) * self.YSize
 
         coordDat = osr.SpatialReference()
         if self.epsgDat:
@@ -285,11 +285,18 @@ def coregister(self,in1,in2):
         trans1 = DS1.GetGeoTransform()
         xsize1 = DS1.RasterXSize
         ysize1 = DS1.RasterYSize
+        epsg1 = self.getProjectionSystem(in1)
 
         DS2 = gdal.Open(in2, gdal.GA_ReadOnly)
         trans2 = DS2.GetGeoTransform()
         xsize2 = DS2.RasterXSize
         ysize2 = DS2.RasterYSize
+        epsg2 = self.getProjectionSystem(in2)
+
+        if epsg1 != epsg2:
+            raise Exception('The current version of geo_autoRIFT assumes the two images are in the same projection, i.e. it cannot handle two different projections; the users are thus recommended to do the tranformation themselves before running geo_autoRIFT.')
+
+
 
         W = np.max([trans1[0],trans2[0]])
         N = np.min([trans1[3],trans2[3]])

diff --git a/geo_autoRIFT/geogrid/src/geogrid.cpp b/geo_autoRIFT/geogrid/src/geogrid.cpp
@@ -1025,8 +1025,8 @@ void geoGrid::geogrid()
 //            {
 //                std::cout << "\n" << lOff+ii << " " << pOff+jj << " " << demLine[jj] << "\n";
 //            }
-            rgind = std::round((rngpix - startingRange) / dr) + 1.;
-            azind = std::round((tline - sensingStart) * prf) + 1.;
+            rgind = std::round((rngpix - startingRange) / dr) + 0.;
+            azind = std::round((tline - sensingStart) * prf) + 0.;
 
 
             //*********************Slant-range vector
@@ -1198,7 +1198,7 @@ void geoGrid::geogrid()
             }
 
 
-            if ((rgind > nPixels)|(rgind < 1)|(azind > nLines)|(azind < 1))
+            if ((rgind > nPixels-1)|(rgind < 1-1)|(azind > nLines-1)|(azind < 1-1))
             {
                 raster1[jj] = nodata_out;
                 raster2[jj] = nodata_out;

diff --git a/geo_autoRIFT/geogrid/src/geogridOptical.cpp b/geo_autoRIFT/geogrid/src/geogridOptical.cpp
@@ -921,8 +921,8 @@ void geoGridOptical::geogridOptical()
                 targutm0[pp] = xyzs[pp];
             }
 
-            xind = std::round((targutm0[0] - startingX) / XSize) + 1.;
-            yind = std::round((targutm0[1] - startingY) / YSize) + 1.;
+            xind = std::round((targutm0[0] - startingX) / XSize) + 0.;
+            yind = std::round((targutm0[1] - startingY) / YSize) + 0.;
 
 
 
@@ -1000,7 +1000,7 @@ void geoGridOptical::geogridOptical()
             }
 
 
-            if ((xind > nPixels)|(xind < 1)|(yind > nLines)|(yind < 1))
+            if ((xind > nPixels-1)|(xind < 1-1)|(yind > nLines-1)|(yind < 1-1))
             {
                 raster1[jj] = nodata_out;
                 raster2[jj] = nodata_out;