6

I have images taken from plane, it's longitude, latitude, elevation, omega, phi, kappa. I suppose that earth is flat (I don't need any DEM on that stage) and I want to calculate how image is nailed to the ground. I want to find center and 4 corners of nailed image to the ground.

Could you advice some simple formula for calculation? Or free tool that has such calculations just to check my own implementation results. I tried to use GRASS but without any luck, it is too complicated :(

Thanks in advance!

UPDATE: Images were taken from around 3.2k meters without physical GCPs. The precision that we got is around 20 meters with camera angles (Omega, Phi) around 2 grad. Do you think it is possible to get better precision? We are using C# and here is our implementation:

    /// <summary>
    /// 
    /// </summary>
    /// <param name="omega">camera angle, in grads</param>
    /// <param name="phi">camera angle, in grads</param>
    /// <param name="kappa">camera angle, in grads</param>
    /// <param name="longitude">camera position</param>
    /// <param name="latitude">camera position</param>
    /// <param name="altitude">camera position</param>
    /// <param name="altOverGround">average flight height, it is defined on flight planning stage</param>
    /// <param name="imageLength">it is defined on flight planning stage, calculated from altOverGround and focal lenght, it is in meters</param>
    /// <param name="imageWidth">it is defined on flight planning stage, calculated from altOverGround and focal lenght, it is in meters</param>
    /// <returns></returns>
    public static DbGeography CalculateProjection(double omega, double phi, double kappa, double longitude, double latitude, double altitude,
        double altOverGround, double imageLength, double imageWidth)
    {
        // coordinates of the camera
        DbGeography cameraCoordinates = CreatePoint(longitude, latitude);

        // converting angles from grad to radian
        omega = ConvertFromGradToRadian(omega);
        phi = ConvertFromGradToRadian(phi);
        kappa = ConvertFromGradToRadian(kappa);

        // get image center (projected to groud) as offset to camera coordinates
        // result of imageCenter is in image coordinate system
        Point imageCenter = GetImageCenter(omega, phi, kappa, altitude);

        // calculating real image length for this height
        var realImageLength = imageCenter.Z * imageLength / altOverGround;
        var realImageWidth = imageCenter.Z * imageWidth / altOverGround;

        // corners offset in image coordinates system
        var xOfsfet = realImageLength / 2;
        var yOffset = realImageWidth / 2;

        // corners list
        List<Point> imageCorners = new List<Point>() { 
            new Point(xOfsfet, yOffset, imageCenter.Z),
            new Point(-xOfsfet, yOffset, imageCenter.Z),
            new Point(-xOfsfet, -yOffset, imageCenter.Z),
            new Point(xOfsfet, -yOffset, imageCenter.Z)
        };

        List<DbGeography> projectionCorners = new List<DbGeography>();
        foreach (var point in imageCorners)
        {
            // get image corner projected to the ground
            Point projectedPoint = GetProjectionCoordinates(omega, phi, kappa, altitude, point);

            // add projected point to list of projected corners
            projectionCorners.Add(CreatePointWithOffset(cameraCoordinates, projectedPoint.Y, projectedPoint.X));
        }

        // add first point to complete the polygon
        projectionCorners.Add(projectionCorners[0]);

        return CreatePolygonFromListOfPoints(projectionCorners);            
    }

    private static Point GetProjectionCoordinates(double omega, double phi, double kappa, double altitude, Point imageCorner)
    {
        // converting point coordinates from image's coordinate system to normal coordinate system
        Point imagePoint = ConvertImageCoordToNormal(omega, phi, kappa, imageCorner);

        // get coordinates of the projected point
        return GetProjectionCoordinates(altitude, imagePoint);
    }

    private static Point GetProjectionCoordinates(double altitude, Point imagePoint)
    {
        var x = altitude * imagePoint.X / imagePoint.Z;
        var y = altitude * imagePoint.Y / imagePoint.Z;
        return new Point(x, y, altitude);
    }

    private static Point ConvertImageCoordToNormal(double omega, double phi, double kappa, Point imageCorner)
    {
        double x1, x2, y1, y2, z1, z2, x, y, z;
        kappa = Math.PI / 2 - kappa;
        x1 = imageCorner.X * Math.Cos(kappa) + imageCorner.Y * Math.Sin(kappa);
        y1 = -imageCorner.X * Math.Sin(kappa) + imageCorner.Y * Math.Cos(kappa);
        z1 = imageCorner.Z;

        x2 = x1 * Math.Cos(phi) - z1 * Math.Sin(phi);
        y2 = y1;
        z2 = x1 * Math.Sin(phi) + z1 * Math.Cos(phi);

        x = x2;
        y = y2 * Math.Cos(omega) + z2 * Math.Sin(omega);
        z = -y2 * Math.Sin(omega) + z2 * Math.Cos(omega);
        return new Point(x, y, z);
    }

    /// <summary>
    /// Calculates altitude of the image center
    /// </summary>
    /// <param name="omega"></param>
    /// <param name="phi"></param>
    /// <param name="kappa"></param>
    /// <param name="altitude"></param>
    /// <returns>Altitude of image center</returns>
    private static Point GetImageCenter(double omega, double phi, double kappa, double altitude)
    {
        // calculating coordinates of image center
        var xCenter = altitude * Math.Tan(phi);
        var yCenter = altitude * Math.Tan(omega);
        var imageAltitude = altitude / (Math.Cos(omega) * Math.Cos(phi));

        // recalculating center with kappa value
        var alpha = Math.Atan(yCenter / xCenter) + kappa;
        var aa = Math.Sqrt(xCenter * xCenter + yCenter * yCenter);

        double imageCenterX = Math.Cos(alpha) * aa;
        double imageCenterY = Math.Sin(alpha) * aa;
        return new Point(imageCenterX, imageCenterY, imageAltitude);
    }

    private static double ConvertFromGradToRadian(double omega)
    {
        omega = Math.PI * omega / 200;
        return omega;
    }

and Point class:

class Point
{
    public double X { get; set; }
    public double Y { get; set; }
    public double Z { get; set; }
    public Point(double x, double y, double z)
    {
        X = x;
        Y = y;
        Z = z;
    }
}
1

Here is Easiest and free way to georeference an image file using QGIS and here are more links

Georeferencing an image in QGis is a very simple task. Here is how:

Run QGis and create a new, blank project. You don't need to save it. You should already have some other georeferenced material, the best is aerial photos, but also shapes could do. Open the folder where those files are, and drop them onto the blank project. QGis will load and show them. In the QGis window, click on Plugins - Georeferencer - Georeferencer. The georeferencing tool will start up and show a window. Click on File - Open raster, and browse to the image file, in my case a jpg file. The plugin loads the file on its window. Click on the Add point tool from the toolbar (it's the icon with three red dots), then click on a known point on the map. It needs to be a point you can recognize on your other map, like the edge of a house, or a special natural feature, etc. A popup appears, asking you to input the coordinates. Click on the "From map" button. The popup temporarily closes, and you are reverted to the original map. Navigate with the keyboard and mouse wheel and locate the point on the map corresponding to the point you clicked on in the previous point. Click there, and you are returned to the popup window. A red dot is placed on the point you have just georeferenced. You repeat points 5 and 6 on several different points. You should choose points on different areas on the map, and not close together. The further away, the better. At least three are necessary. Once you finished, click on File - Start Georeferencing. A popup appears. You should fill at least the "Output raster" file. It will generate a geofererenced tiff file. You can also play with the "Transformation type" box. For me, for example, the Helmert transform worked best. Click on OK and the plugin will process your file. Now you should close the window (you can save your points in a separate file, in case you need to repeat the process later, which I recommend doing). Finally, you should check your map. A quick and easy way is to use the main window, which already contains your reference map. Now, drag the tiff file you generated on the main window: the map will be added. Double click on its layed on the left pane, and move the Global Transparency slider to somewhere aroud 50%. Now you can compare the original map and the new one, to see if it matches. If is does not, you can go back to the previous points and adjust some parameters. Your tiff file is now georeferenced, using the same coordinate system as the map you used as reference. It is now ready to be imported in mapserver.

  • I don't have GCP so first of all I need nail image to the ground by calculating image projection with Omega, Phi, Kappa. I am still looking forward for precise formula or tool that can do that. After that I will have 4 corner points of image. And yes I will use gdal(qgisis using gdal for georeferencing) to georeference my image, that should not be a problem. – Serhiy Zhdynyak Dec 26 '13 at 11:35
  • I saw that document and much more other articles and documents. Based on that we implement our own solution that tries to reproject image to the ground, but the result are not precise, and I am looking for tools or some kind of final formula to check my results. So I will know that my implementation works ok. Anyway thanks for helping! – Serhiy Zhdynyak Dec 26 '13 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.