NITF IGEOLO bounds interpretation

[warmerdam]

Hi Frank,

Just noticed there is a bug in the NITF IGEOLO handling. According to the
NITF spec, the  IGEOLO field in the NITF header contains the four corner
coordinates in the following order:
(0,0),(0,MaxCol),(MaxRow,MaxCol),(MaxRow,0) where MaxCol = NCOLS - 1 and
MaxRow = NROWS - 1. Note: MaxCol, MaxRow represents the topleft corner of
the bottomright pixel, not the bottomright corner of the bottomright pixel.
Consequently, I've made the following changes to nitfdataset.cpp:

1. On read, changed the following code:
---------------------------------
       poDS->adfGeoTransform[1] =
           (psImage->dfLRX - psImage->dfULX) / poDS->nRasterXSize;
       poDS->adfGeoTransform[2] = 0.0;
       poDS->adfGeoTransform[3] = psImage->dfULY;
       poDS->adfGeoTransform[4] = 0.0;
       poDS->adfGeoTransform[5] =
           (psImage->dfLRY - psImage->dfULY) / poDS->nRasterYSize;
---------------------------------

to

---------------------------------
       poDS->adfGeoTransform[1] =
           (psImage->dfLRX - psImage->dfULX) / (poDS->nRasterXSize - 1);
       poDS->adfGeoTransform[2] = 0.0;
       poDS->adfGeoTransform[3] = psImage->dfULY;
       poDS->adfGeoTransform[4] = 0.0;
       poDS->adfGeoTransform[5] =
           (psImage->dfLRY - psImage->dfULY) / (poDS->nRasterYSize - 1);
---------------------------------

2. On write, changed the following code:

---------------------------------
   dfURX = dfULX + padfGeoTransform[1] * nRasterXSize;
   dfURY = dfULY + padfGeoTransform[4] * nRasterXSize;
   dfLRX = dfULX + padfGeoTransform[1] * nRasterXSize
       + padfGeoTransform[2] * nRasterYSize;
   dfLRY = dfULY + padfGeoTransform[4] * nRasterXSize
       + padfGeoTransform[5] * nRasterYSize;
   dfLLX = dfULX + padfGeoTransform[2] * nRasterYSize;
   dfLLY = dfULY + padfGeoTransform[5] * nRasterYSize;
---------------------------------

to

---------------------------------
   dfURX = dfULX + padfGeoTransform[1] * (nRasterXSize - 1);
   dfURY = dfULY + padfGeoTransform[4] * (nRasterXSize - 1);
   dfLRX = dfULX + padfGeoTransform[1] * (nRasterXSize - 1)
       + padfGeoTransform[2] * (nRasterYSize - 1);
   dfLRY = dfULY + padfGeoTransform[4] * (nRasterXSize - 1)
       + padfGeoTransform[5] * (nRasterYSize - 1);
   dfLLX = dfULX + padfGeoTransform[2] * (nRasterYSize - 1);
   dfLLY = dfULY + padfGeoTransform[5] * (nRasterYSize - 1);
---------------------------------

David Carter

[warmerdam]

David,

I have reviewed the NITF 2.1 specification and I am still uncertain of
the correct interpretation.  The NITF 2.1 spec says specifically:

"""The locations of the four corners of the (significant) image data shall 
be given in image coordinate order: (0,0), (0, MaxCol), (MaxRow,
MaxCol), (MaxRow, 0). MaxCol and MaxRow shall be determined 
from the values contained, respectively, in NCOLS and NROWS 
as MaxCol = NCOLS -1 and MaxRow = NROWS -1."""

So the locations are referred to as being corners of the significant image
data.  The listed pixel locations for the corners seems to suggest these
locations are given for ordering purposes, and it isn't clear to me that
the "corners of the significant imagery" at found at the top left corner
of the pixels indicated.  Can you point me to some additional material 
that would support your interpretation? 

I would add that the corner coordinates are sometimes initialized
from the RPF location header information.  I have reviewed 
MIL-STD-2411-2 to see what the details are for it's corner coordinates
but I was unable to find any useful detail text in that document.  But
if we change the interpretation of IGEOLO we may have to do it
differently to avoid inadvertant re-interpretation of the RPF location
information.  

I tried looking at the IGEOLO and RPF location information for some
datasets, but due to the lack of precision of the IGEOLO data I was
unable to reach any meaningful conclusions other than that the RPF
corner coordinates *seem* to refer to the outside edge of the imagery
data based on the tiles abutting without overlap based on that 
interpretation. 

To ensure we don't lose track of this I have created a bugzilla entry for
it.  I would suggest we append further information there for posterities
sake. 

Best regards,

Frank

[warmerdam]

2. Changed the following lines in NITFDataset::Open() (nitfdataset.cpp):

-----------------------------------
       poDS->bGotGeoTransform = TRUE;
       poDS->adfGeoTransform[0] = psImage->dfULX;
       poDS->adfGeoTransform[1] =
           (psImage->dfLRX - psImage->dfULX) / poDS->nRasterXSize;
       poDS->adfGeoTransform[2] = 0.0;
       poDS->adfGeoTransform[3] = psImage->dfULY;
       poDS->adfGeoTransform[4] = 0.0;
       poDS->adfGeoTransform[5] =
           (psImage->dfLRY - psImage->dfULY) / poDS->nRasterYSize;
-----------------------------------

to:

-----------------------------------
       poDS->bGotGeoTransform = TRUE;
       poDS->adfGeoTransform[1] = (psImage->dfLRX - psImage->dfULX) /
(poDS->nRasterXSize - 1);
       poDS->adfGeoTransform[5] = (psImage->dfLRY - psImage->dfULY) /
(poDS->nRasterYSize - 1);
       poDS->adfGeoTransform[0] = psImage->dfULX - 0.5 * poDS->adfGeoTransform[1];
       poDS->adfGeoTransform[2] = 0.0;
       poDS->adfGeoTransform[3] = psImage->dfULY - 0.5 * poDS->adfGeoTransform[5];
       poDS->adfGeoTransform[4] = 0.0;
-----------------------------------

3.  Changed the following lines in NITFDataset::SetGeoTransform()
(nitfdataset.cpp):

-----------------------------------
   dfULX = padfGeoTransform[0];
   dfULY = padfGeoTransform[3];
   dfURX = dfULX + padfGeoTransform[1] * nRasterXSize;
   dfURY = dfULY + padfGeoTransform[4] * nRasterXSize;
   dfLRX = dfULX + padfGeoTransform[1] * nRasterXSize
       + padfGeoTransform[2] * nRasterYSize;
   dfLRY = dfULY + padfGeoTransform[4] * nRasterXSize
       + padfGeoTransform[5] * nRasterYSize;
   dfLLX = dfULX + padfGeoTransform[2] * nRasterYSize;
   dfLLY = dfULY + padfGeoTransform[5] * nRasterYSize;
-----------------------------------

to:

-----------------------------------
   dfULX = padfGeoTransform[0] + 0.5 * padfGeoTransform[1] + 0.5 *
padfGeoTransform[2];
   dfULY = padfGeoTransform[3] + 0.5 * padfGeoTransform[4] + 0.5 *
padfGeoTransform[5];
   dfURX = dfULX + padfGeoTransform[1] * (nRasterXSize - 1);
   dfURY = dfULY + padfGeoTransform[4] * (nRasterXSize - 1);
   dfLRX = dfULX + padfGeoTransform[1] * (nRasterXSize - 1)
                       + padfGeoTransform[2] * (nRasterYSize - 1);
   dfLRY = dfULY + padfGeoTransform[4] * (nRasterXSize - 1)
                       + padfGeoTransform[5] * (nRasterYSize - 1);
   dfLLX = dfULX + padfGeoTransform[2] * (nRasterYSize - 1);
   dfLLY = dfULY + padfGeoTransform[5] * (nRasterYSize - 1);
-----------------------------------

David Carter
Software Developer
ER Mapper

[warmerdam]

Hi Frank,

I received a response to my question about the IGEOLO field in NITF
subheader (see below).  As stated below, the NITF standard does not address
my question, but the Compendium of Controlled Extensions, Image Chip (ICHIP)
specification does.  It specifies the center of the pixel as the
registration point.  So I'm inclined to use the pixel center, but I guess
it's not that important because the IGEOLO field is only approximate and not
intended for analytical purposes.

David Carter
Software Developer
ER Mapper

[warmerdam]

Hi Frank,

I noticed the bugzilla bug report did not include the contents of the email
(see attached) I received from Bridget Durham (JITC NITF Lab) in response to
my question.  To quote Bridget:

"...the NITF standard itself does not address this question, but the
Compendium of Controlled Extensions, Image Chip (ICHIP) specification does.
It specifies the center of the pixel as this point."

This was the main reason I thought it would be better to use the pixel
centre.  I also reasoned that if you use the pixel centre, you'll be at most
1/2 pixel out if some other application associates the IGEOLO coordinates
with the pixel corners.

[warmerdam]

I am now in the situation where I am *generally* willing to treat the IGEOLO
as the center of the corner pixels. 

However, I am pretty confident that the RPF CoverageSectionSubheader is
in terms of outer pixel edges, and it is also transmitted via the NITFImage
dfULX/etc fields.  So the IGEOLO changes cannot safely be applied as shown.  
 
I think the solution would be to still treat the psImage bounds as they are 
now (as outer edges), and to apply adjustments in the code that reads and
writes from IGEOLO.

[warmerdam]

David's IGEOLO changes have now been incorporated (in altered form) into 
GDAL.