Iammuratc / satellitepy

Fine-grained object detection in satellite images
MIT License
2 stars 0 forks source link

Gaofen dataset merge subimages #109

Open UnknownHobbyist opened 1 year ago

UnknownHobbyist commented 1 year ago

Merging the subimages back into the original images and updating bounding boxes accordingly.

UnknownHobbyist commented 1 year ago

A more general approach of doing the same task without the stitching class from OpenCV. We are using the latitute and longitute (standard wgs84 format) from the no dublicates files). We are only depending on the utm and pyproj modules which is necessary to obtain the x and y coordinate from latitute and longitute. Also see the last line which is needed to sort the existing no_dublicates files.

See below.
UnknownHobbyist commented 1 year ago

One thing that should be noted is that there are flaws in the no_dublicates file. In the following example you can see that the spatial_resolution (probably in relation to the wrong wkt_epsg_4326 values) is way to low which naturally results in incorrect variables in the script. The spatial resolution for these images must be manually examined.

{
        "sequence_id": "593db4f4-f3ee-4acb-a11a-a24032182d58",
        "base_images": [
            {
                "id": "5a601627-01a9-4016-adc0-25542652c5f4",
                "image_id": "879470e2-745d-4f78-bdce-4f81a956b008",
                "image_path": "/home/murat/Projects/airplane_recognition/data/Gaofen/test/images/1.tif",
                "image_width": 1024,
                "image_height": 1024,
                "spatial_resolution_dy": -5.36441802978516e-06,
                "spatial_resolution_dx": 5.36441802978516e-06,
                "orig_projection": "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]",
                "orig_wkt": "POLYGON ((103.9361572265625 30.5680847167969, 103.941650390625 30.5680847167969, 103.941650390625 30.5625915527344, 103.9361572265625 30.5625915527344, 103.9361572265625 30.5680847167969))",
                "wkt_epsg_4326": "POLYGON ((103.9361572265625 30.5680847167969, 103.941650390625 30.5680847167969, 103.941650390625 30.5625915527344, 103.9361572265625 30.5625915527344, 103.9361572265625 30.5680847167969))",
                "ground_truth": [
                    {
                        "pixel_position": "POLYGON ((909 959, 1031 989, 1056 889, 933 859, 909 959))",
                        "orig_projection": "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]",
                        "orig_wkt": "POLYGON ((103.9410334825516 30.56294023990634, 103.9416879415512 30.56277930736544, 103.941822052002 30.56331574916842, 103.9411622285843 30.56347668170931, 103.9410334825516 30.56294023990634))",
                        "wkt_epsg_4326": "POLYGON ((103.9410334825516 30.56294023990634, 103.9416879415512 30.56277930736544, 103.941822052002 30.56331574916842, 103.9411622285843 30.56347668170931, 103.9410334825516 30.56294023990634))",
                        "class": "A330"
                    },
                    {
                        "pixel_position": "POLYGON ((800 586, 765 657, 825 687, 860 616, 800 586))",
                        "orig_projection": "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]",
                        "orig_wkt": "POLYGON ((103.9404487609863 30.56494116783145, 103.9402610063553 30.56456029415133, 103.9405828714371 30.56439936161044, 103.9407706260681 30.56478023529055, 103.9404487609863 30.56494116783145))",
                        "wkt_epsg_4326": "POLYGON ((103.9404487609863 30.56494116783145, 103.9402610063553 30.56456029415133, 103.9405828714371 30.56439936161044, 103.9407706260681 30.56478023529055, 103.9404487609863 30.56494116783145))",
                        "class": "A321"
                    }
                ]
            }
        ]
    }
Iammuratc commented 1 year ago

Great job! Could you please save the stitched images and the labels on 2tb-0?

UnknownHobbyist commented 1 year ago

You can find the images under /mnt/2tb-0/satellitepy/data/Gaofen/merged/images along with the no_dublicates_merged.json which provides the updated bounding boxes.

UnknownHobbyist commented 1 year ago

Make the no_dublicates_merged.json in satellitepy format.

UnknownHobbyist commented 1 year ago

Manually fixed the wrong spatial resolution problem. Works with the newest script.

from satellitepy.utils.path_utils import get_project_folder
from pathlib import Path
from satellitepy.data.labels import *
from pyproj import Transformer

import json
import numpy as np
import cv2 as cv
import math
import utm

def main():
    imgs = []
    img_files = []

    folder_path = Path("/mnt/2tb-0/satellitepy/data/Gaofen")
    sequences_file = folder_path / "sequences.json"

    no_duplicates_file_train = folder_path / "train/recognition/no_duplicates_sorted.json"
    no_duplicates_file_val = folder_path / "val/recognition/no_duplicates_sorted.json"
    no_duplicates_file_test = folder_path / "test/recognition/no_duplicates_sorted.json"

    f = open(str(sequences_file), "r")
    sequences = json.loads(f.read())
    f.close()

    f1 = open(str(no_duplicates_file_train), "r")
    no_duplicates_train = json.loads(f1.read())
    f1.close()

    f2 = open(str(no_duplicates_file_val), "r")
    no_duplicates_val = json.loads(f2.read())
    f2.close()

    f3 = open(str(no_duplicates_file_test), "r")
    no_duplicates_test = json.loads(f3.read())
    f3.close()

    create_folder(get_project_folder() / "in_folder/Gaofen")
    create_folder(get_project_folder() / "in_folder/Gaofen/train")
    create_folder(get_project_folder() / "in_folder/Gaofen/train/img")
    create_folder(get_project_folder() / "in_folder/Gaofen/train/labels")

    create_folder(get_project_folder() / "in_folder/Gaofen/test")
    create_folder(get_project_folder() / "in_folder/Gaofen/test/img")
    create_folder(get_project_folder() / "in_folder/Gaofen/test/labels")

    create_folder(get_project_folder() / "in_folder/Gaofen/val")
    create_folder(get_project_folder() / "in_folder/Gaofen/val/img")
    create_folder(get_project_folder() / "in_folder/Gaofen/val/labels")

    transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857")

    my_list = []

    counter = {
        'train' : 0,
        'val' : 0,
        'test' : 0
    }

    for i in range(0, len(sequences)):

        labels = init_satellitepy_label()
        available_tasks = ['obboxes', 'fine-class']
        all_tasks = get_all_satellitepy_keys()
        not_available_tasks = [task for task in all_tasks if not task in available_tasks or available_tasks.remove(task)]

        ids = [int(str(folder_path / x['relative_path']).split("/")[-1].split(".")[0])-1 for x in sequences[i]['images']]
        img_files = [str(folder_path / x['relative_path']) for x in sequences[i]['images']]

        if sequences[i]['images'][0]['relative_path'].startswith('train'):
            no_duplicates = no_duplicates_train
            category = "train"
        elif sequences[i]['images'][0]['relative_path'].startswith('val'):
            no_duplicates = no_duplicates_val
            category = "val"
        else:
            no_duplicates = no_duplicates_test
            category = "test"

        min_x = 100000000
        max_x = 0
        min_y = 100000000
        max_y = 0

        xy_corner = []
        flag = False

        for id in ids:
            poly = no_duplicates[id]['base_images'][0]['wkt_epsg_4326']
            poly = poly.replace("POLYGON ((", "").replace("))", "").split(", ")

            wgs_bbox = [[float(y) for y in x.split(" ")] for x in poly]

            s_res_y = 1 / -no_duplicates[id]['base_images'][0]['spatial_resolution_dy']
            s_res_x = 1 / no_duplicates[id]['base_images'][0]['spatial_resolution_dx']

            if not no_duplicates[id]['base_images'][0]['spatial_resolution_dx'] == 0.81:
                xy_bbox = [[int(transformer.transform(x[1], x[0])[0] * s_res_x), int(transformer.transform(x[1], x[0])[1] * s_res_y)] for x in wgs_bbox]
            else:
                xy_bbox = [[int(utm.from_latlon(x[1], x[0])[0] * s_res_x), int(utm.from_latlon(x[1], x[0])[1] * s_res_y)] for x in wgs_bbox]

            xy_min = [min([x[0] for x in xy_bbox]), min([x[1] for x in xy_bbox])]

            min_x = min(min_x, abs(xy_min[0]))
            max_x = max(max_x, abs(xy_min[0]))
            min_y = min(min_y, abs(xy_min[1]))
            max_y = max(max_y, abs(xy_min[1]))

            xy_corner.append(xy_min)

        if flag:
            counter[category] += 1
            continue

        max_x = max_x - min_x + 1024
        max_y = max_y - min_y + 1024

        stitch = np.zeros((max_y, max_x, 3)).astype('uint8')

        for j in range(len(img_files)):

            img = cv.imread(img_files[j])
            x = xy_corner[j][0] - min_x
            y = xy_corner[j][1] - np.sign(xy_corner[j][1]) * min_y

            if xy_corner[j][1] < 0:
                stitch[-y:-y+1024, x:x+1024, 0:3] = img
            else:
                stitch[max_y - y - 1024:max_y - y, x:x+1024, 0:3] = img

            ground_truth = no_duplicates[ids[j]]['base_images'][0]['ground_truth']

            for label in ground_truth:
                bbox = [[int(y) for y in x.split(" ")] for x in label['pixel_position'].replace("POLYGON ((", "").replace("))", "").split(", ")[:-1]]

                if xy_corner[j][1] < 0:
                    bbox = [[int(x + corner[0]), int(-y + corner[1])] for corner in bbox]
                else: 
                    bbox = [[int(x + corner[0]), int((max_y - y - 1024) + corner[1])] for corner in bbox]

                if not is_dublicate(labels['obboxes'], bbox):
                    labels['obboxes'].append(bbox)
                    labels['fine-class'].append(label['class'])

                fill_none_to_empty_keys(labels, not_available_tasks)

        out_file = get_project_folder() / f"in_folder/Gaofen/{category}/img/img_{counter[category]}.tif"
        cv.imwrite(str(out_file), stitch)
        print(f"Image {i} successfully merged into ", out_file)

        f1 = open(str(get_project_folder() / f"in_folder/Gaofen/{category}/labels/{counter[category]}.json"), "w")
        f1.write(json.dumps(labels, indent=4))
        f1.close()

        counter[category] += 1

def is_dublicate(bboxes, searchle):

    for bbox in bboxes:
        error = 0

        for coord1, coord2 in zip(bbox, searchle):

            error += int(math.sqrt((coord1[0] - coord2[0])**2 + (coord1[1] - coord2[1])**2))

        if error < 30:
            return True

    return False

def create_folder(folder):

    Path(folder).mkdir(parents=True, exist_ok=True)

if __name__ == '__main__':
    main()

    # following function provides the sorted json files for no dublicates
    # no_duplicates_val = sorted(no_duplicates_val, key = lambda item: int(Path(item['base_images'][0]['image_path']).stem.split(".")[0]))