rasoberry-noaa/rectify.py

#!/usr/bin/python3

from multiprocessing import Pool, cpu_count
import sys
import re
import numpy
from math import atan,sin,cos,sqrt,tan,acos,ceil
from PIL import Image

EARTH_RADIUS = 6371.0
SAT_HEIGHT = 830.0
SAT_ORBIT_RADIUS = EARTH_RADIUS + SAT_HEIGHT
SWATH_KM = 2800.0
THETA_C = SWATH_KM / EARTH_RADIUS

# Note: theta_s is the satellite viewing angle, theta_c is the angle between the projection of the satellite on the
# Earth's surface and the point the satellite is looking at, measured at the center of the Earth

# Compute the satellite angle of view given the center angle


def theta_s(theta_c):
    return atan(EARTH_RADIUS * sin(theta_c)/(SAT_HEIGHT+EARTH_RADIUS*(1-cos(theta_c))))

# Compute the inverse of the function above


def theta_c(theta_s):
    delta_sqrt = sqrt(EARTH_RADIUS**2 + tan(theta_s)**2 *
                      (EARTH_RADIUS**2-SAT_ORBIT_RADIUS**2))
    return acos((tan(theta_s)**2*SAT_ORBIT_RADIUS+delta_sqrt)/(EARTH_RADIUS*(tan(theta_s)**2+1)))

# The nightmare fuel that is the correction factor function.
# It is the reciprocal of d/d(theta_c) of theta_s(theta_c) a.k.a.
# the derivative of the inverse of theta_s(theta_c)


def correction_factor(theta_c):
    norm_factor = EARTH_RADIUS/SAT_HEIGHT
    tan_derivative_recip = (
        1+(EARTH_RADIUS*sin(theta_c)/(SAT_HEIGHT+EARTH_RADIUS*(1-cos(theta_c))))**2)
    arg_derivative_recip = (SAT_HEIGHT+EARTH_RADIUS*(1-cos(theta_c)))**2/(EARTH_RADIUS*cos(
        theta_c)*(SAT_HEIGHT+EARTH_RADIUS*(1-cos(theta_c)))-EARTH_RADIUS**2*sin(theta_c)**2)

    return norm_factor * tan_derivative_recip * arg_derivative_recip

# Radians position given the absolute x pixel position, assuming that the sensor samples the Earth
# surface with a constant angular step


def theta_center(img_size, x):
    ts = theta_s(THETA_C/2.0) * (abs(x-img_size/2.0) / (img_size/2.0))
    return theta_c(ts)

# Worker thread


def wthread(rectified_width, corr, endrow, startrow):
    # Make temporary working img to push pixels onto
    working_img = Image.new(img.mode, (rectified_width, img.size[1]))
    rectified_pixels = working_img.load()

    for row in range(startrow, endrow):
        # First pass: stretch from the center towards the right side of the image
        start_px = orig_pixels[img.size[0]/2, row]
        cur_col = int(rectified_width/2)
        target_col = cur_col

        for col in range(int(img.size[0]/2), img.size[0]):
            target_col += corr[col]
            end_px = orig_pixels[col, row]
            delta = int(target_col) - cur_col

            # Linearly interpolate
            for i in range(delta):
                # For night passes of Meteor the image is just gray level and
                # start_px and end_px being an int instead of a tuple
                if type(start_px) != int:
                    interp_r = int((start_px[0]*(delta-i) + end_px[0]*i) / delta)
                    interp_g = int((start_px[1]*(delta-i) + end_px[1]*i) / delta)
                    interp_b = int((start_px[2]*(delta-i) + end_px[2]*i) / delta)
                    rectified_pixels[cur_col,row] = (interp_r, interp_g, interp_b)
                else:
                    interp = int((start_px*(delta-i) + end_px*i) / delta)
                    rectified_pixels[cur_col,row] = interp
                cur_col += 1

            start_px = end_px

        # First pass: stretch from the center towards the left side of the image
        start_px = orig_pixels[img.size[0]/2, row]
        cur_col = int(rectified_width/2)
        target_col = cur_col

        for col in range(int(img.size[0]/2)-1, -1, -1):
            target_col -= corr[col]
            end_px = orig_pixels[col, row]
            delta = cur_col - int(target_col)

            # Linearly interpolate
            for i in range(delta):
                # For night passes of Meteor the image is just gray level and
                # start_px and end_px being an int instead of a tuple
                if type(start_px) != int:
                    interp_r = int((start_px[0]*(delta-i) + end_px[0]*i) / delta)
                    interp_g = int((start_px[1]*(delta-i) + end_px[1]*i) / delta)
                    interp_b = int((start_px[2]*(delta-i) + end_px[2]*i) / delta)
                    rectified_pixels[cur_col,row] = (interp_r, interp_g, interp_b)
                else:
                    interp = int((start_px*(delta-i) + end_px*i) / delta)
                    rectified_pixels[cur_col,row] = interp
                cur_col -= 1

            start_px = end_px

    # Crop the portion we worked on
    slice = working_img.crop(box=(0, startrow, rectified_width, endrow))
    # Convert to a numpy array so STUPID !#$&ING PICKLE WILL WORK
    out = numpy.array(slice)
    # Make dict of important values, return that.
    return {"offs": startrow, "offe": endrow, "pixels": out}


if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Usage: {} <input file>".format(sys.argv[0]))
        sys.exit(1)

    out_fname = re.sub("\..*$", "-rectified", sys.argv[1])

    img = Image.open(sys.argv[1])
    print("Opened {}x{} image".format(img.size[0], img.size[1]))

    # Precompute the correction factors
    corr = []
    for i in range(img.size[0]):
        corr.append(correction_factor(theta_center(img.size[0], i)))

    # Estimate the width of the rectified image
    rectified_width = ceil(sum(corr))

    # Make new image
    rectified_img = Image.new(img.mode, (rectified_width, img.size[1]))

    # Get the pixel 2d arrays from the source image
    orig_pixels = img.load()

    # Callback function to modify the new image
    def modimage(data):
        if data:
            # Write slice to the new image in the right place
            rectified_img.paste(Image.fromarray(
                data["pixels"]), box=(0, data["offs"]))

    # Number of workers to be spawned - Probably best to not overdo this...
    numworkers = cpu_count()
    # Estimate the number of rows per worker
    wrows = ceil(img.size[1]/numworkers)
    # Initialize some starting data
    startrow = 0
    endrow = wrows
    # Make out process pool
    p = Pool(processes=numworkers)

    # Let's have a pool party! Only wnum workers are invited, though.
    for wnum in range(numworkers):
        # Make the workers with appropriate arguments, pass callback method to actually write data.
        p.apply_async(wthread, (rectified_width, corr,
                                endrow, startrow), callback=modimage)
        # Aparrently ++ doesn't work?
        wnum = wnum+1
        # Beginning of next worker is the end of this one
        startrow = wrows*wnum
        # End of the worker is the specified number of rows past the beginning
        endrow = startrow + wrows
        # Show how many processes we're making!
        print("Spawning process ", wnum)
    # Pool's closed, boys
    p.close()
    # It's a dead pool now
    p.join()

    rectified_img.save(out_fname + ".png", "PNG")