[gdal-dev] NITF Chip

[Even Rouault]

Le lundi 10 juin 2013 22:43:52, Baker, Anthony W a écrit :
> One possible answer seems to be to use ICHIP_OP_**_** and  ICHIP_FI_**_**
> 
> Using the _OP_ metadata I form a row-wise 2x2 projection matrix with _12 -
> _11 as one row and _21 - _11 as the other.  With _FI_ I form a column-wise
> 2x2 projection matrix with _12 - _11 and _21 - _11 as the columns.  
> Applying these matrices and the two translations in the right order seems
> to give ok results.  If I combine this transform with the geoTransform
> (matrix multi and trans) I seem to get ok results.  The results are
> different by a little from ArcMap's raster inclusion but close.  Code is
> included below to show the transforms that I have done.  Does this seem
> right?
> 
> double fi11r = strtod(poDataset->GetMetadataItem("ICHIP_FI_ROW_11", NULL),
> NULL); double fi11c = strtod(poDataset->GetMetadataItem("ICHIP_FI_COL_11",
> NULL), NULL); double fi12r =
> strtod(poDataset->GetMetadataItem("ICHIP_FI_ROW_12", NULL), NULL); double
> fi12c = strtod(poDataset->GetMetadataItem("ICHIP_FI_COL_12", NULL), NULL);
> double fi21r = strtod(poDataset->GetMetadataItem("ICHIP_FI_ROW_21", NULL),
> NULL); double fi21c = strtod(poDataset->GetMetadataItem("ICHIP_FI_COL_21",
> NULL), NULL); double op11r =
> strtod(poDataset->GetMetadataItem("ICHIP_OP_ROW_11", NULL), NULL); double
> op11c = strtod(poDataset->GetMetadataItem("ICHIP_OP_COL_11", NULL), NULL);
> double op12r = strtod(poDataset->GetMetadataItem("ICHIP_OP_ROW_12", NULL),
> NULL); double op12c = strtod(poDataset->GetMetadataItem("ICHIP_OP_COL_12",
> NULL), NULL); double op21r =
> strtod(poDataset->GetMetadataItem("ICHIP_OP_ROW_21", NULL), NULL); double
> op21c = strtod(poDataset->GetMetadataItem("ICHIP_OP_COL_21", NULL), NULL);
> 
> //Column matrix
> double pr11, pr12, pr22, pr21;
> double npr1, npr2;
> npr1 = sqrt((fi12c - fi11c)*(fi12c - fi11c) + (fi12r - fi11r)*(fi12r -
> fi11r)); npr2 = sqrt((fi21c - fi11c)*(fi21c - fi11c) + (fi21r -
> fi11r)*(fi21r - fi11r)); pr11 = (fi12c - fi11c)/npr1;
> pr12 = (fi21c - fi11c)/npr2;
> pr21 = (fi12r - fi11r)/npr1;
> pr22 = (fi21r - fi11r)/npr2;
> 
> //Row matrix
> double br11, br12, br22, br21;
> double nbr1, nbr2;
> nbr1 = sqrt((op12c - op11c)*(op12c - op11c) + (op12r - op11r)*(op12r -
> op11r)); nbr2 = sqrt((op21c - op11c)*(op21c - op11c) + (op21r -
> op11r)*(op21r - op11r)); br11 = (op12c - op11c)/nbr1;
> br21 = (op21c - op11c)/nbr2;
> br12 = (op12r - op11r)/nbr1;
> br22 = (op21r - op11r)/nbr2;
> 
> //Transform
> double m11, m12, m22, m21;
> m11 = pr11*br11 + pr12*br21;
> m21 = pr21*br11 + pr22*br21;
> m12 = pr11*br12 + pr12*br22;
> m22 = pr21*br12 + pr22*br22;
> 
> //Compute the translation factor
> double a1, a2;
> a1 = -(m11*op11c + m12*op11r) + fi11c;
> a2 = -(m21*op11c + m22*op11r) + fi11r;
> 
> ST[0] = a1;
> ST[3] = a2;
> ST[1] = m11;
> ST[2] = m12;
> ST[4] = m21;
> ST[5] = m22;

I've not followed the math logic of the above formulas, but if I were you I 
would try defining 4 GDAL GCPs with dfGCPPixel = OP_COL, dfGCPLine = OP_ROW, 
dfGCPX = FI_COL, dfGCPY = FI_ROW and then I would use GDALGCPsToGeoTransform() 
(that will do a least squares error regression) to deduce the best 
geotransform matrix.


-- 
Geospatial professional services
http://even.rouault.free.fr/services.html