Revision: 6625
          http://sourceforge.net/p/jump-pilot/code/6625
Author:   michaudm
Date:     2020-11-19 08:02:57 +0000 (Thu, 19 Nov 2020)
Log Message:
-----------
Fix more deeply and document the half pixel shift problem

Modified Paths:
--------------
    
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoImage.java
    
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoReferencedRaster.java

Modified: 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoImage.java
===================================================================
--- 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoImage.java   
    2020-11-15 17:54:16 UTC (rev 6624)
+++ 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoImage.java   
    2020-11-19 08:02:57 UTC (rev 6625)
@@ -43,6 +43,7 @@
 import java.awt.image.renderable.ParameterBlock;
 import java.lang.reflect.InvocationTargetException;
 
+import javax.media.jai.InterpolationNearest;
 import javax.media.jai.JAI;
 import javax.media.jai.RenderedOp;
 import javax.media.jai.operator.AffineDescriptor;
@@ -147,7 +148,7 @@
 //          System.out.println("GI: NO SCALE CACHE");
 
           // First, scale the original image
-          double scaleX = scale * gtr.getDblModelUnitsPerRasterUnit_X();
+          double scaleX = scale * 
Math.abs(gtr.getDblModelUnitsPerRasterUnit_X());
           double scaleY = scale * 
Math.abs(gtr.getDblModelUnitsPerRasterUnit_Y());
 
           // calculate predicted dimensions
@@ -209,7 +210,11 @@
             // Interpolation interp = Interpolation
             // .getInstance(Interpolation.INTERP_BICUBIC);
             // pb.add(interp); // add interpolation method
-            img = JAI.create("scale", pb, hints);
+            //img = JAI.create("scale", pb, hints);
+            // Prefer affine to scale as affine uses double parameters and 
preserve scale signus
+            // while scale uses float parameters and forbid negative scaleY
+            AffineTransform tr = new AffineTransform(scaleX_toUse, 0, 0, 
scaleY_toUse, 0, 0);
+            img = JAI.create("affine", scale_src_img, tr, new 
InterpolationNearest());
           } else {
             pb.add(scaleX_toUse);
             pb.add(scaleY_toUse);
@@ -270,9 +275,9 @@
       raster_cropX = Math.min((float)raster_cropX, (float)img.getWidth());
       raster_cropY = Math.min((float)raster_cropY, (float)img.getHeight());
       raster_cropW = Math
-          .min((float)raster_cropW, (float)img.getWidth() - /*(int)*/ 
raster_cropX);
-      raster_cropH = Math.min((float)raster_cropH, (float)img.getHeight()
-          - /*(int)*/ raster_cropY);
+          .min((float)raster_cropW, (float)(img.getWidth() - raster_cropX));
+      raster_cropH = Math.min((float)raster_cropH, (float)(img.getHeight()
+          - raster_cropY));
 
       pb = new ParameterBlock();
       pb.addSource(img);
@@ -282,8 +287,11 @@
       //System.out.println("croph " + (float)raster_cropH + " " + 
(img.getHeight() - /*(int)*/ raster_cropY));
       pb.add((float)raster_cropX);
       pb.add((float)raster_cropY);
-      pb.add((float)raster_cropW);
-      pb.add((float)raster_cropH);
+      // conversions of cropX/Y and width/height from double to float may
+      // create an envelope slightly outside the real envelope of the image
+      // removing ulp to the width/height solve the problem
+      pb.add((float)raster_cropW-Math.ulp((float)raster_cropW));
+      pb.add((float)raster_cropH-Math.ulp((float)raster_cropH));
       img = JAI.create("crop", pb, null);
 
       // move the image to the model coordinates

Modified: 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoReferencedRaster.java
===================================================================
--- 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoReferencedRaster.java
    2020-11-15 17:54:16 UTC (rev 6624)
+++ 
core/trunk/src/com/vividsolutions/jump/workbench/imagery/geoimg/GeoReferencedRaster.java
    2020-11-19 08:02:57 UTC (rev 6625)
@@ -41,7 +41,6 @@
 import java.io.IOException;
 import java.io.InputStream;
 import java.net.URI;
-import java.net.URISyntaxException;
 import java.util.Arrays;
 import java.util.List;
 
@@ -77,6 +76,10 @@
   @Deprecated
   Coordinate coorModel_tiepointLT;
 
+  // AreaOrPoint=AREA (default) means that image coordinates
+  // refer to the upper left corner of the upper left angle.
+  // AreaOrPoint=POINT means that image coordinates refer to
+  // the center of the upper left pixel
   public enum AreaOrPoint {AREA,POINT}
   AreaOrPoint areaOrPoint = AreaOrPoint.AREA;
 
@@ -91,6 +94,8 @@
   private Coordinate modelULPixelCenter;
 
   private double dblModelUnitsPerRasterUnit_X;
+  // [michaudm 2020-11-16] signed y-scale
+  // to be able to handle north-south oriented or south-north oriented images
   private double dblModelUnitsPerRasterUnit_Y;
 
   /**
@@ -120,7 +125,7 @@
       re = new ReferencedImageException("probably no tiff image: "
           + e.getMessage());
     } catch (IOException e) {
-      re = new ReferencedImageException("problem acessing tiff image: "
+      re = new ReferencedImageException("problem accessing tiff image: "
           + e.getMessage());
     } finally {
       // clean up
@@ -155,7 +160,7 @@
       tags[4] = fieldModelGeoTransform.getAsDouble(3);
       // y-ordinate of the center of the upper left pixel
       tags[5] = fieldModelGeoTransform.getAsDouble(7);
-      Logger.debug("gtiff trans: " + Arrays.toString(tags));
+      Logger.debug("gtiff transform: " + Arrays.toString(tags));
       setEnvelope(tags);
     }
     // use the tiepoints as defined
@@ -166,25 +171,20 @@
       // Map the modeltiepoints from raster to model space
 
       // Read the tiepoints
-      double offset = (areaOrPoint == AreaOrPoint.AREA) ? 0.5 : 0.0;
-      // in areaOrPoint=AREA model, reference for the TiePoint is the UL
-      // corner of the pixel
-      // To set the transformation parameters we use the center of the
-      // UL coordinate
+      // imageCoord has not the same meaning in AREA and in POINT mode,
+      // but here, we don't mind, imageCoord may represent any point in
+      // the image, important thing is that imageCoord and modelCoord
+      // represent the same point in the image and in the model
+      // coordinate system.
       Coordinate imageCoord = new Coordinate(
-              fieldModelTiePoints.getAsDouble(0)-offset,
-              fieldModelTiePoints.getAsDouble(1)-offset);
+              fieldModelTiePoints.getAsDouble(0),
+              fieldModelTiePoints.getAsDouble(1));
       Coordinate modelCoord = new Coordinate(
               fieldModelTiePoints.getAsDouble(3),
               fieldModelTiePoints.getAsDouble(4));
 
-      //setCoorRasterTiff_tiepointLT(new Coordinate(
-      //    fieldModelTiePoints.getAsDouble(0),
-      //    fieldModelTiePoints.getAsDouble(1), 0));
-      //setCoorModel_tiepointLT(new Coordinate(
-      //    fieldModelTiePoints.getAsDouble(3),
-      //    fieldModelTiePoints.getAsDouble(4), 0));
-      //setEnvelope();
+      Logger.debug("gtiff tiepoints found : " + 
Arrays.toString(fieldModelTiePoints.getAsDoubles()));
+
       // Find the ModelPixelScale field
       XTIFFField fieldModelPixelScale = dir
           .getField(XTIFF.TIFFTAG_GEO_PIXEL_SCALE);
@@ -195,8 +195,6 @@
             + "\n" + MSG_GENERAL);
       }
 
-      Logger.debug("gtiff tiepoints found : " + 
Arrays.toString(fieldModelTiePoints.getAsDoubles()));
-
       dblModelUnitsPerRasterUnit_X = fieldModelPixelScale.getAsDouble(0);
       dblModelUnitsPerRasterUnit_Y = fieldModelPixelScale.getAsDouble(1);
       //https://trac.osgeo.org/gdal/ticket/4977
@@ -203,14 +201,20 @@
       if (!honourNegativeScaleY) dblModelUnitsPerRasterUnit_Y = - 
Math.abs(dblModelUnitsPerRasterUnit_Y);
       Logger.debug("gtiff scale : scalex=" + dblModelUnitsPerRasterUnit_X + ", 
scaley=" + dblModelUnitsPerRasterUnit_Y);
 
-      double tx = 
modelCoord.x-getDblModelUnitsPerRasterUnit_X()*(imageCoord.x);
-      double ty = 
modelCoord.y-getDblModelUnitsPerRasterUnit_Y()*(imageCoord.y);
+      // To compute the translation parameters of the transformation, we need
+      // to know how a point is converted from image to model coordinate, we
+      // don't mind where this point is.
+      double tx = modelCoord.x - dblModelUnitsPerRasterUnit_X * imageCoord.x;
+      double ty = modelCoord.y - dblModelUnitsPerRasterUnit_Y * imageCoord.y;
 
-      // Compute the model coordinate for the center of UL and LR pixels
-      rasterULPixelCenter = new Coordinate(0.0,0.0);
+      // Now, we want to know the model coordinate of the point precisely
+      // located at the center of the upper left pixel.
+      // Coordinates of this point is not the same in AREA and in POINT modes
+      rasterULPixelCenter = (areaOrPoint == AreaOrPoint.AREA) ?
+              new Coordinate(0.5,0.5) : new Coordinate(0.0,0.0);
       modelULPixelCenter = new Coordinate(
-              getDblModelUnitsPerRasterUnit_X()*rasterULPixelCenter.x + tx,
-              getDblModelUnitsPerRasterUnit_Y()*rasterULPixelCenter.y + ty);
+              getDblModelUnitsPerRasterUnit_X() * rasterULPixelCenter.x + tx,
+              getDblModelUnitsPerRasterUnit_Y() * rasterULPixelCenter.y + ty);
 
       setEnvelope();
 
@@ -294,7 +298,7 @@
         tags[i] = Double.parseDouble(line);
       }
 
-      Logger.debug("wf: " + Arrays.toString(tags));
+      Logger.debug("worldfile: " + Arrays.toString(tags));
 
       setEnvelope(tags);
 
@@ -355,24 +359,17 @@
 
     // set up a default envelope
     double[] tags = new double[6];
-    tags[0] = 1; // pixel size in x
-                 // direction
+    tags[0] = 1; // pixel size in x direction
     tags[1] = 0; // rotation about y-axis
     tags[2] = 0; // rotation about x-axis
-    tags[3] = -1;// pixel size in the
-                 // y-direction
-    tags[4] = 0; // x-coordinate of the
-                 // center of the upper
-                 // left pixel
-    tags[5] = 0; // y-coordinate of the
-                 // center of the upper
-                 // left pixel
+    tags[3] = -1;// pixel size in the y-direction
+    tags[4] = 0; // x-coordinate of the center of the upper left pixel
+    tags[5] = 0; // y-coordinate of the center of the upper left pixel
     setEnvelope(tags);
   }
 
   private void setEnvelope(double[] tags) {
-    //setCoorRasterTiff_tiepointLT(new Coordinate(-0.5, -0.5));
-    //setCoorModel_tiepointLT(new Coordinate(0, 0));
+
     AffineTransform transform = new AffineTransform(tags);
 
     //We should honour negative scale y, but gdal created plenty of
@@ -386,11 +383,14 @@
     dblModelUnitsPerRasterUnit_X = transform.getScaleX();
     dblModelUnitsPerRasterUnit_Y = transform.getScaleY();
 
-    Point2D rasterLT = new Point2D.Double(src.getMinX(), src.getMinY());
+    // To compute the envelope in AreaOrPoint.AREA mode, we need to
+    // know that upper left pixel center is at 0.5, 0.5, not 0, 0
+    double offset = (areaOrPoint == AreaOrPoint.AREA) ? 0.5 : 0.0;
+    Point2D rasterLT = new Point2D.Double(src.getMinX()+offset, 
src.getMinY()+offset);
     Point2D modelLT = new Point2D.Double();
     transform.transform(rasterLT, modelLT);
 
-    rasterULPixelCenter = new Coordinate(src.getMinX(), src.getMinY());
+    rasterULPixelCenter = new Coordinate(rasterLT.getX(), rasterLT.getY());
     modelULPixelCenter = new Coordinate(modelLT.getX(), modelLT.getY());
 
     setEnvelope();
@@ -397,22 +397,23 @@
   }
 
   void setEnvelope() {
-    double offset = (areaOrPoint == AreaOrPoint.AREA) ? 0.5 : 0.0;
-    // Get the image coordinate of the bottom left corner of the bottom left 
pixel
-    // from the image coordinate of the center of the bottom left pixel
-    Coordinate coorRaster_imageLB = new Coordinate(
-            rasterULPixelCenter.x-offset,src.getHeight()-offset);
-    // Get the image coordinate of the top right corner of the top right pixel
-    // from the image coordinate of the center of the top right pixel
-    Coordinate coorRaster_imageRT = new Coordinate(
-            src.getWidth()-offset, -offset);
-    // Transform the bottom left and the top right corner to model coordinate
-    // to get the envelope of the image
-    Coordinate coorModel_imageLB = rasterToModelSpace(coorRaster_imageLB);
-    Coordinate coorModel_imageRT = rasterToModelSpace(coorRaster_imageRT);
 
-    envModel_image = new Envelope(coorModel_imageLB, coorModel_imageRT);
+    // Image coordinate of the upper left corner of the envelope
+    double ulx = rasterULPixelCenter.x-0.5;
+    double uly = rasterULPixelCenter.y-0.5;
 
+    // Bottom left coordinate of the envelope
+    Coordinate imageEnvelopeBL = new Coordinate(ulx,uly+src.getHeight());
+
+    // Top right coordinate of the envelope
+    Coordinate imageEnvelopeTR = new Coordinate(ulx+src.getWidth(), uly);
+
+    // Transform envelope corners to the model coordinate system
+    Coordinate modelEnvelopeBL = rasterToModelSpace(imageEnvelopeBL);
+    Coordinate modelEnvelopeTR = rasterToModelSpace(imageEnvelopeTR);
+
+    envModel_image = new Envelope(modelEnvelopeBL, modelEnvelopeTR);
+
     // backup original envelope
     envModel_image_backup = envModel_image;
   }
@@ -426,13 +427,14 @@
    */
   private Coordinate rasterToModelSpace(Coordinate coorRaster) {
     Coordinate coorModel = new Coordinate();
-
     coorModel.x = modelULPixelCenter.x
         + (coorRaster.x - rasterULPixelCenter.x)
         * dblModelUnitsPerRasterUnit_X;
+    System.out.println("" + modelULPixelCenter.x + " + (" + coorRaster.x + "-" 
+ rasterULPixelCenter.x + ")*" + dblModelUnitsPerRasterUnit_X + " = " + 
coorModel.x);
     coorModel.y = modelULPixelCenter.y
         + (coorRaster.y - rasterULPixelCenter.y)
         * dblModelUnitsPerRasterUnit_Y;
+    System.out.println("" + modelULPixelCenter.y + " + (" + coorRaster.y + "-" 
+ rasterULPixelCenter.y + ")*" + dblModelUnitsPerRasterUnit_Y + " = " + 
coorModel.y);
     coorModel.z = 0;
 
     return coorModel;



_______________________________________________
Jump-pilot-devel mailing list
Jump-pilot-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/jump-pilot-devel

Reply via email to