Foreword
The mosaic line is actually the seam in the field of cv. When doing remote sensing image mosaic, in order to make the transition of the splicing natural, it is necessary for the line at the splicing to present a winding state. Some needs also require the mosaic line to avoid houses, take the terrain line and so on. Mosaic line algorithm, in ENVI the auto generate seamlines can be very intuitive to see the generation. However, the mosaic lines generated in ENVI are very straight, directly passing through some obvious curve textures, which makes the mosaic image leave obvious mosaic marks. The realization of this algorithm is mainly to solve the problem of straight splicing line.
Introduction to the algorithm
The main points mainly include several parts. Generally, mosaic is carried out after orthographic correction. After orthographic correction, the image will be rotated to a certain extent, and the black area generated is the area filled with 0 value. This part of the area does not need to be processed, so we need to extract the effective area when processing; The method of generating paths based on image texture features is mainly obtained by using the super-pixel segmentation method in the previous article; the general idea of mosaic line mosaic is to generate lines first, then cut out the irregular edges of a single image, and finally splice them according to the warp method given by gdal.
Introduction to packages used
| whl name | Version | whl name | version | whl name | version | whl name | version |
|---|---|---|---|---|---|---|---|
| os | basic whl | cv2 | 4.6.0.66 | networkx | 2.8.6 | pathlib | basic whl |
| math | basic whl | ogr | 3.4.3 | itertools | basic whl | geopandas | 0.11.1 |
| shapefile | basic whl | osr | 3.4.3 | shapely | 2.3.1 | skimage | 0.19.3 |
| tarfile | basic whl | gdal | 3.4.3 | shutil | basic whl | numpy | 1.23.3 |
| time | basic whl | gdalconst | 3.4.3 | topojson | 1.5 | python | 3.10.4 |
Some of the names of the packages are not quite right, but most of them are correct, mainly scikit-image and networkx gdal and opencv. Some of the interesting points are that the shapely package needs to be written when it is installed
|
|
Full code
The code is very long, below according to the different function, take the main function as the mentality, carries on disassembles one by one.
Master function
|
|
In order to facilitate debugging and checking, the process files are generated one by one here, but in fact, gdal has a corresponding memory method that virtul raster can directly put the array into memory, which will be rewritten as a memory method later when the project is close to completion.
Grid binarization
|
|
Here mainly aims at the multi-band grid, has carried on the binarization to extract the grid effective scope, uses is the cv package.
After binarization, only 0 and 1 are left in the grid to facilitate grid to vector assignment.
Binary grid to vector
|
|
The binarization grid also mainly uses the osr library. Those who are familiar with gdal should know that gdal, OSR, ogr and other geographic data processing packages have been osgeo gradually incorporated into the management since a certain year, using the open source foundation to generate the power of sustainable development. The grid value of the previous binarization is used to polygon facilitate the subsequent extraction of the SHP block with a value of 1.
Select a polygon block based on the attribute values
|
|
Based on the ogr package, the vector block with DN value equal to 1 is selected through the key function SetAttributeFilter, and a new SHP is created to be saved
Find the processing area of the mosaic line
This is mainly to use the effective area SHP to find the intersection. After the above functions are used to process the two images to be inlaid, the intersection of the two images is obtained based on the geopandas package.
|
|
Of course, after getting the intersection, part of the boundary of the vector is jagged, because the process of converting the grid to the vector does not simplify the jaggies. The reduced function is given here.
|
|
Clips the imported original grid with the valid area
|
|
This is mainly based on the warp function to clip the original grid, obtain the overlapping area of the two images and output it as TIF, which is convenient for the subsequent use of super-resolution segmentation.
Superpixel segmentation
Because SLIC quickshift it takes a long time to process large-resolution images, I use downsampling to reduce the resolution of images, and then use low-resolution images to segment them. The test speed is increased from 500 seconds to 90 seconds. The resampling function is as follows:
|
|
Then there is the main segmentation function.
|
|
The SLIC segmentation algorithm is used here. The parameter settings are basically conventional, and there is not much attempt. The reason is that this is simply to obtain the network path. In fact, how complex it is is useless. If it is too square, it is not as good as Voronoi that, so it is set in the middle.
Segmentation result processing
Convert the split grid to a vector line. Difficulties have been encountered here, because the gdal grid to vector can only be converted from the grid to the face, but not directly to the line, so it is necessary to convert the grid to the face first and then to the line:
|
|
The code is also a conventional way of writing, each step of the operation will be output to a file to facilitate the inspection and analysis of the problem. After the segmentation result line is output, a simplification is performed to remove the sawtooth problem of the gdal raster to vector. A new method is used here to simplify the function:
|
|
As a result of the vector line after division, each block is a separate object, and we need to merge all objects into one object:
|
|
The test found that there are still some inexplicable scattered lines in the result of merging, and sometimes there are various problems, so in order to be on the safe side, reverse thinking here is to intersect the dividing line with the dividing area to get the overall line layer, which contains only one object:
|
|
Gets the start and end points of the tessellation line
The idea here is to take the diagonal points of the SHP blocks in the grid intersection area as the starting and ending points, and use vector packages such as geopandas to easily obtain the corner coordinates:
|
|
Some people 
here may ask why the corner is not at the edge of the image. In fact, the overall range of the image is like this: 
only because of the setting of arcmap’s ignore background value, the background of the image display is ignored.
Optimal Path Algorithm
In the article posted some time ago, there is actually an introduction to this, which is mainly based on the networkx package. There are astar three options for the algorithm, dijkstra, and bellman-ford. Here is a star:
|
|
Vector line clipping face
Unexpectedly, this is still very difficult. At present, there should be very few bloggers to write this introduction. Here is a brief introduction of ideas. The purpose is to use the generated optimal path polyline to clip the range polygon of the grid. First, the line is used as a buffer, and the buffer value needs to be very small. Then, the generated buffer polygon surface grid range SHP, which is almost a line, is used to calculate the difference, and the split SHP can be obtained:
|
|
Although the correlation function is also written in the previous text, the number of water words is given here.
After processing the above content, in fact, the two faces to be split are still in the same object. We need to “explode” the integrated object, then sort it according to the size of the area, and select the area with the largest area to output, so that we can get the real mosaic line clipping area of the grid.
|
|
Inlay
The result of the clipping of the tessellation line is obtained, and then the formal tessellation is performed:
|
|
Here I also studied, it seems that gdal can not directly write the option in a warp in one step, and then cut the direct mosaic, can only be done step by step, basically is the warp function, gdal is really powerful.
Here, due to the better quality of the two images, there is no clear dividing line.
Post-processing work
Mainly refers to the work of the pyramid, because it is the final work, so the previous generation process did not build the pyramid, I only build the pyramid in the mosaic image, mainly for the convenience of viewing, this code is not difficult.
|
|
Sum up
This idea is also running-in for a long time. At the beginning, I thought about using cv seamless stitch, but I found that it would destroy the real value of the image. I could only use physical methods to generate mosaic lines. The task requirement is to avoid the house. There is no idea here. What I thought before was to use the mask function of super-pixel segmentation to make a mask for the house area. Based on the NDVI and OTSU threshold method, the building area is extracted as the non-selected area for segmentation, but it is found that the optimal path will still follow all areas, but based on arcmap, it can normally follow the non-building area. The reason may be that the underlying algorithm of ArcGIS is different from that of gdal, but I don’t have much ability and energy to study the underlying problems, so I can only divide them according to the architectural texture as much as possible at present.
The main code is to change the file output of vector and raster. It is time-consuming to read and write all the time. In practical application, it should be rewritten as a virtual dataset unified processing package for vector. However, there is also a problem involved here, such as how to operate the vectorization of virtual raster, such as the conversion of raster to surface. Existing packages are generated based on local raster files. If they are placed in a virtual raster, will they be stored in an array? This is an unknown problem that needs to be verified.
Although geopandas is powerful, there are too many dependent packages, including Fiona, but in fact, Fiona and other packages can also implement all the above algorithms about vectors. It remains to be verified whether the vector technology stack can be replaced by ogr or Fiona.
The generation network path mentioned voronoi above is actually an alternative to SLIC, because the generation speed of Voronoi diagram is theoretically far superior to SLIC or quickshift, but the triangulation generated by Voronoi diagram is too square, which will lead to the loss of some texture routes. If we add the manual editing of mosaic lines in the later stage, that is, the semi-automatic generation of mosaic lines, in fact, the idea of Voronoi diagram should be desirable. I have not studied the generation of Voronoi diagram due to time constraints, and we need to discuss and understand it in the future.
Appendix · General Code
|
|