A batch texture mapping method for irregular top surface of 3D model

Abstract

This application discloses a batch texture mapping method for irregular top surfaces of 3D models, including: S1. Generating matching texture images based on GIS features; S2. Calculating the rotation angle of the minimum bounding box of the GIS features; S3. Adding the rotation angle of the minimum bounding box to the GIS features; S4. Creating texture mapping CGA rules in CityEngine; S5. Applying texture mapping to the irregular top surface of the 3D model and correcting it. This application enables batch texture mapping of the top surface of 3D models. Based on GIS vector and raster data, the generated texture images can strictly maintain consistency with the shape of each model's top surface, thus ensuring the accuracy of the texture files; at the same time, by calculating and adjusting the texture rotation angle, the mapping accuracy is further guaranteed. As an automated batch processing solution, compared with traditional manual texture mapping, this application significantly saves labor costs and improves work efficiency.

Description

Technical Field

[0001] This application belongs to the field of 3D modeling technology of geographic information systems, specifically involving a batch texture mapping method for irregular top surfaces of 3D models. Background Technology

[0002] With the rapid development of smart city, digital twin, and Geographic Information System (GIS) technologies, large-scale three-dimensional (3D) modeling is increasingly widely used in urban planning, emergency management, and virtual reality. In high-precision 3D scenes, the quality of the texture mapping on the model surface directly determines the realism and visualization effect of the scene. Among them, the top surface of the 3D model (such as building roofs, terrain surfaces, etc.) is one of the visual focal points, and the clarity, fit, and geometric consistency of its texture are crucial to the overall modeling quality.

[0003] However, in large-scale 3D modeling of GIS scenes, texture mapping on the top surface of the model faces significant complexity and technical challenges. Unlike planar textures of the ground, the top surface of a 3D model often exhibits high geometric complexity, containing numerous non-rectangular irregular polygons with varying orientations and rotation angles. This diversity of geometric forms makes it difficult to directly apply general-purpose automatic projection algorithms, easily leading to problems such as texture stretching, misalignment, obvious seams, or inverted orientation, making it difficult to meet the requirements of high-precision modeling.

[0004] To address the aforementioned issues, current methods for processing top surface textures primarily rely on manual operations within 3D modeling software (such as 3ds Max, Maya, Blender, CityEngine, etc.). Since top surfaces often involve irregular polygons, manual texturing remains the most effective and quality-assured method in current industry practice. Specifically, technicians need to manually perform operations such as UV unwrapping, texture coordinate adjustment, image cropping, and rotation matching to ensure that the texture image matches the shape of the model's top surface.

[0005] While manual texturing ensures accuracy for individual models to some extent, it has significant shortcomings and limitations in large-scale batch modeling applications. Firstly, it is inefficient and difficult to scale. Handling massive amounts of 3D models at the city level requires extremely long processing times, making it difficult to meet project delivery requirements. Secondly, it incurs high labor costs. Texturing work heavily relies on the experience of professional modelers and involves repetitive work, leading to a significant increase in project labor budgets. Finally, it suffers from poor consistency and maintainability. Manual operation is susceptible to subjective factors, making it difficult to ensure consistent texturing standards across batches of models. Furthermore, the maintenance cost of re-texturing when the underlying GIS data is updated is extremely high.

[0006] In summary, overcoming the problems of low efficiency and high cost of traditional manual texturing, while solving the technical bottlenecks of shape matching and angle correction in automatic texturing of irregular polygon top surfaces, and achieving high-precision, batch-based automatic texturing based on GIS data, has become an urgent technical problem to be solved in the field of GIS 3D modeling. Summary of the Invention

[0007] This application provides a batch texture mapping method for irregular top surfaces of 3D models, in order to solve the technical problems of low efficiency, poor consistency and maintainability, and difficulty in scaling up when performing texture mapping on irregular shapes of the top surface of 3D models in large-scene 3D modeling technology of existing geographic information systems due to reliance on manual operation.

[0008] To address the aforementioned technical problems, this application adopts the following technical solution: a batch texture mapping method for irregular top surfaces of 3D models, comprising:

[0009] S1. Generate matching texture images based on GIS features;

[0010] S2. Calculate the rotation angle of the smallest rectangular bounding box of a GIS feature;

[0011] S3. Add the rotation angle of the minimum rectangular bounding box to the GIS feature;

[0012] S4. Create a texture mapping CGA rule in CityEngine;

[0013] S5. Apply texture mapping and correct the irregular top surface of the 3D model.

[0014] Furthermore, the method in step S1 includes:

[0015] S11. Prepare GIS vector data shp_file and raster data tif_file that are fully corresponding to the coordinate system. Configure the texture image output folder output_dir, the field name id_field used as feature identifier in shp_file, and the invalid raster value no_data_value in tif_file. The default value of no_data_value is 0.

[0016] S12. Import the Arcpy and PIL.Image libraries, set arcpy.env.overwriteOutput = True, and check if output_dir, shp_file, tif_file and id_field exist in turn. If they do not exist, throw an exception.

[0017] S13. Use the arcpy.da.SearchCursor() function to read the id_field and feature geometry object field "SHAPE@" from shp_file, create a cursor and convert it into a feature list, traverse the feature list, and obtain the feature ID and feature_shape of each feature.

[0018] S14. Call the arcpy.management.MakeRasterLayer() function to register tif_file as a temporary in-memory raster map raster_layer, and use the arcpy.management.Clip() function to clip raster_layer according to feature_shape to obtain a temporary clipped raster file temp_raster;

[0019] S15. Use the arcpy.RasterToNumPyArray() function to convert temp_raster to a NumPy array raster_array. Based on the dimension of raster_array, classify and process it to generate an RGBA format PNG texture image with an alpha channel. Output a prompt message if the image is invalid. Delete temp_raster after traversing all elements.

[0020] Furthermore, the method in step S2 includes:

[0021] S21. Define the cross product function cross(), which takes the coordinates of the starting point O, point A, and point B as input and returns the cross product value of the two-dimensional vectors pointing from O to A and from O to B.

[0022] S22. Define a convex hull function convex_hull(), take a point set points as input, sort the point set in ascending order of x-coordinate, and sort in ascending order of y-coordinate if the x-coordinates are the same, determine the direction of the points by the cross product function, construct the lower convex hull and the upper convex hull respectively, remove duplicate endpoints and concatenate to obtain the convex hull vertex sequence and return it.

[0023] S23. Define the minimum bounding rectangle function min_bounding_rectangle(), take the convex hull point set hull_points as input, traverse each edge of the convex hull, calculate the direction angle of the current edge and rotate hull_points in the reverse direction until the current edge is parallel to the X-axis, calculate the area of ​​the axis-aligned bounding box of the rotated point set, filter out the bounding box with the smallest area, restore its corner points to the original coordinate system, and return the vertex sequence of the minimum bounding rectangle;

[0024] S24. Define a calculation function calc_mbr_and_angles(), take a polygon point set polygon_points as input, call the convex hull function and the minimum bounding box function in sequence to obtain the minimum bounding rectangle mbr of the feature, calculate the direction angles angle1 and angle2 of two adjacent sides of mbr with due east as 0 degrees and the angles increasing counterclockwise as the rule, and return mbr, angle1 and angle2 in dictionary form.

[0025] Furthermore, the specific rules for determining the direction of a point using the cross product function in step S22 are as follows:

[0026] If the cross product is less than 0 after a new point is added, it indicates that the path turns to the right, and the vertex is removed from the stack; if the cross product is equal to 0, it indicates that the points are collinear, and the vertex is removed from the stack; if the cross product is greater than 0, it indicates that the path turns to the left, and the new point is retained and added to the stack.

[0027] Furthermore, the specific steps for calculating the directional angles of two adjacent sides of mbr in step S24 are as follows:

[0028] Extract the coordinates of adjacent vertices of the MBR and calculate the vectors edge1 and edge2 of two adjacent edges. Calculate the radian values ​​of edge1 and edge2 using the math.atan2() function, and then convert the radian values ​​to angle values ​​using the math.degrees() function to obtain angle1 and angle2.

[0029] Furthermore, the method in step S3 includes:

[0030] S31. Prepare the GIS vector data shp_file to which the rotation angle will be added, and set the initial name of the field to be added to 'rangle';

[0031] S32. Import the Arcpy library, check if the field 'rangle' exists in shp_file. If it does not exist, add a 'rangle' field of type FLOAT using the arcpy.AddField_management() function. If it exists, update the field name to 'rotate_angle' and add a field of the same type. Denote the newly added field as output_field.

[0032] S33. Read the "FID", "SHAPE@" and output_field fields of shp_file through the arcpy.da.UpdateCursor() function, create an update cursor and traverse it to obtain the feature number fid and polygon object polygon, and extract the set of vertices polygon_points of polygon;

[0033] S34. Call the calc_mbr_and_angles() function to get angle1 and angle2, select the angle with the smaller absolute value as the feature rotation angle rangle, assign float(rangle) to the output_field field, and update the feature attribute record through cursor.updateRow(row). After traversing all cursors, the rotation angle is added.

[0034] Furthermore, the method in step S4 includes:

[0035] S41. Create a CGA rule file topSkin.cga, use the attr keyword to declare the attribute Id as the texture index, initialize the value to 0, and define the img() function to dynamically return the file path and name of the corresponding texture image based on the Id;

[0036] S42. Use the `attr` keyword to declare the attribute `rangle` as the texture rotation angle, with an initial value of 0; declare the boolean attribute `isRotate` to control the secondary rotation of the texture, with an initial value of false, and explicitly define its optional values ​​as true and false using the `@Enum` annotation; define the function `calcRotateAngle()`, which returns 180 - `rangle` if `isRotate` is true, and returns - `rangle` if `isRotate` is false.

[0037] S43. Use the attr keyword to declare the boolean properties xFlip, yFlip, and zFlip that control the mirror flipping of the X, Y, and Z axes. The default configuration is xFlip=true, yFlip=false, and zFlip=false.

[0038] S44. Create the TopSkin core rule for top surface texture mapping. Perform coordinate system alignment, rotation, and mirroring in sequence to generate a quadrilateral consistent with the original shape bounding box of the feature. Set the texture projection layer and coordinate system, and call the Skin rule to perform the mapping operation. The Skin rule first checks whether the texture image corresponding to img(Id) exists through the fileExists() function. If it exists, it performs the texture (img(Id)) mapping operation. Otherwise, it performs the NIL operation to delete the current shape.

[0039] Furthermore, in step S44,

[0040] The specific CGA operation commands for performing coordinate system alignment, rotation, and mirroring are as follows:

[0041] First, alignScopeToAxes(y) is executed to align the current shape's range coordinate system to the global Y-axis. Then, alignScopeToGeometryBBox(xz) is executed to align the coordinate system to the XZ plane of the target geometry's bounding box. Next, rotateScope(0, calcRotateAngle, 0) is executed to rotate the coordinate system around the Y-axis by the angle calcRotateAngle. Finally, mirrorScope(xFlip, yFlip, zFlip) is executed to complete the coordinate system mirroring process.

[0042] The specific CGA operation commands for generating quadrilaterals and setting texture projection are as follows:

[0043] The primitiveQuad() operation generates a quadrilateral that matches the bounding box for the original shape; the setupProjection(0, scope.xz, scope.sx, scope.sz) and projectUV(0) operations are executed in sequence to complete the texture layer setup and projection coordinate system creation.

[0044] Furthermore, the method in step S5 includes:

[0045] S51. Import the GIS vector data shp_file with Id and rangle attributes into CityEngine, generate a shape layer and verify the integrity of the attributes, generate a 3D model based on the shape layer, select the top surface shape of all models, apply the TopSkin rule, and generate a rough model with preliminary texture.

[0046] S52. For models where the texture direction does not match the geometry in the rough model, adjust the isRotate property to true in the attribute inspector to rotate the texture by 180 degrees based on the initial rotation angle, thereby completing the precise texture correction and realizing batch texture mapping of irregular top surfaces of 3D models.

[0047] The beneficial effects of this application are as follows: In large-scale 3D modeling of geographic information systems, this application enables batch texturing of the top surface of 3D models. Based on GIS vector and raster data, the generated texture images can strictly maintain consistency with the shape of the top surface of each model (regardless of whether it is regular or not), thereby ensuring the accuracy of the texture files; at the same time, by calculating and adjusting the texture rotation angle, the texturing accuracy is further guaranteed. As an automated batch processing solution, compared with traditional manual texturing, this application significantly saves labor costs and improves work efficiency. Attached Figure Description

[0048] Figure 1This is a flowchart illustrating an embodiment of the batch texture mapping method for irregular top surfaces of 3D models according to this application.

[0049] Figure 2 This is a texture image generated based on vector data, representing an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application.

[0050] Figure 3 This is an example of the minimum bounding rectangle generated based on an irregular polygon, representing an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application.

[0051] Figure 4 This is an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application, showing the addition of a rotation angle field and its assignment effect in the GIS vector data attribute table;

[0052] Figure 5 This is a model rendering with preliminary texture generated according to an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application.

[0053] Figure 6 This is a before-and-after image of manually rotating a texture, representing an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application.

[0054] Figure 7 This is a final texture mapping effect diagram of the top surface shape of a 3D model, according to an embodiment of the batch texture mapping method for irregular top surfaces of 3D models in this application. Detailed Implementation

[0055] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to specific embodiments.

[0056] Numerous specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways than those described herein, and therefore the invention is not limited to the specific embodiments disclosed in the following specification.

[0057] See Figure 1 , Figure 1 This is a flowchart illustrating an embodiment of the batch texture mapping method for irregular top surfaces of 3D models according to this application. The method includes:

[0058] S1. Generate matching texture images based on GIS features.

[0059] Specifically, the method in step S1 includes:

[0060] S11. Prepare GIS vector data and raster data with fully corresponding coordinate systems. The former is denoted as shp_file and the latter as tif_file. Name the texture image output folder as output_dir, the field name used as id in the shp_file file is id_field, and the invalid raster value in the tif_file file is no_data_value, with a default value of 0.

[0061] S12. Import the Arcpy and PIL.Image libraries, and set arcpy's overwrite output to True, i.e., arcpy.env.overwriteOutput = True. Check if the output folder output_dir and the input files shp_file and tif_file exist. If they exist, continue execution; otherwise, throw an exception. Check if the id_field exists in the shp_file. If it exists, continue execution; otherwise, throw an exception.

[0062] S13. Based on the id_field and "SHAPE@" (the geometric object of the feature), use the arcpy.da.SearchCursor() function to read the shp_file and create a cursor, denoted as cursor. Convert the cursor to a list, i.e., rows = list(cursor), and traverse each feature by rows.

[0063] Let the current traversed element be row. Obtain the element number feature_id (i.e. row[0]) and the element geometry feature_shape (i.e. row[1]) based on row.

[0064] S14. Create a temporary layer for the raster data to be clipped. Call the arcpy.management.MakeRasterLayer() function to register the original raster data tif_file as a temporary layer raster_layer in memory.

[0065] Create a temporary cropped raster file, denoted as temp_raster, with an initial value of os.path.join(output_dir, f"temp_{feature_id}.tif"). Use the arcpy.management.Clip() function to crop the raster_layer by feature_shape and save it as temp_raster.

[0066] S15. Convert the cropped raster data into a NumPy array. Use the arcpy.RasterToNumPyArray() function to convert temp_raster to a NumPy array raster_array.

[0067] See Figure 2 Generate PNG images from a raster_array. Use the `shape` attribute of the raster_array to determine the array's dimensions and classify it as a multi-band image, a single-band image, or an invalid image.

[0068] (1) Multi-band images require the raster_array to have a dimension of 3 and at least 3 bands. In this case, the first three bands are extracted as red (red_channel), green (green_channel), and blue (blue_channel) channels respectively. Then, each channel is normalized and converted to uint8 type. Next, np.dstack is used to merge the three normalized channels to generate an array rgb_array. Then, based on the invalid value no_data_value, an alpha channel is set for rgb_array (i.e., the pixels containing the invalid value are set to the alpha channel), and the resulting array is denoted as alpha_channel. Subsequently, the np.dstack() function is used to merge rgb_array and alpha_channel into an RGB array with an alpha channel, denoted as rgba_array. Finally, the Image.fromarray() function is used to generate image data img from rgba_array in 'RGBA' mode, and the img.save() function is used to save it as a png format image. Figure 2 This displays the effect of generating corresponding texture images based on GIS data features.

[0069] (2) For single-band images, the raster_array must have a dimension of 2. Then, the raster_array is normalized, and the resulting data is denoted as normalized. Next, the channel data alpha is created using the method np.where(raster_array == no_data_value, 0, 255).astype(np.uint8). Then, the np.dstack() function is used to merge the normalized and alpha data into an RGB array with an alpha channel, denoted as rgba_array. Finally, the Image.fromarray() function is used to generate image data img from rgba_array in 'RGBA' mode, and the img.save() function is used to save it as a PNG image.

[0070] The core Python code for converting the above two types of images into PNG images is as follows:

[0071] # If the raster is multi-band RGB (e.g., 3-band), process each band.

[0072] if len(raster_array.shape) == 3 and raster_array.shape[0] >= 3:

[0073] # Extract RGB data for each band, ignoring the Alpha channel (if it exists).

[0074] red_channel = raster_array[0, :, :]

[0075] green_channel = raster_array[1, :, :]

[0076] blue_channel = raster_array[2, :, :]

[0077] # Normalize each band of data and convert it to uint8

[0078] def normalize_channel(channel): # Define the internal normalization function

[0079] if channel.dtype == np.uint8: return channel

[0080] min_val = channel.min()

[0081] max_val = channel.max()

[0082] return ((channel - min_val) / (max_val

[0083] - min_val) * 255).astype(np.uint8)

[0084] ####

[0085] red_channel = normalize_channel(red_channel)

[0086] green_channel = normalize_channel(green_channel)

[0087] blue_channel = normalize_channel(blue_channel)

[0088] # Merge RGB bands

[0089] rgb_array = np.dstack((red_channel, green_channel,

[0090] blue_channel))

[0091] # Set the transparency channel (invalid values ​​are transparent), 3 channels RGB

[0092] alpha_channel = np.ones(rgb_array.shape[:2],

[0093] dtype=np.uint8) * 255

[0094] alpha_channel[np.any(raster_array[:3] == no_data_value,

[0095] axis=0)] = 0

[0096] # Merge RGB and Alpha channels

[0097] rgba_array = np.dstack((rgb_array, alpha_channel))

[0098] # Convert an array to an image; 'RGBA' mode includes an alpha channel.

[0099] img = Image.fromarray(rgba_array, 'RGBA')

[0100] # Save as PNG image

[0101] img_path = os.path.join(output_dir, f"{feature_id}.png")

[0102] img.save(img_path)

[0103] elif len(raster_array.shape) == 2: # Single band

[0104] # Normalized data

[0105] min_val, max_val = raster_array.min(), raster_array.max()

[0106] normalized = ((raster_array - min_val) / (max_val

[0107] - min_val) * 255).astype(np.uint8)

[0108] # Create an RGBA image (grayscale + transparent)

[0109] alpha = np.where(raster_array == no_data_value, 0,

[0110] 255).astype(np.uint8)

[0111] rgba_array = np.dstack((normalized, normalized,

[0112] normalized, alpha))

[0113] img = Image.fromarray(rgba_array, 'RGBA')

[0114] img.save(os.path.join(output_dir, f"{feature_id}.png"))

[0115] else:

[0116] print(f"\tFeature ID={feature_id} is not a valid RGB raster!")

[0117] # Clean up temporary files

[0118] if arcpy.Exists(temp_raster):

[0119] arcpy.Delete_management(temp_raster)

[0120] ####

[0121] S2. Calculate the rotation angle of the smallest rectangular bounding box of a GIS feature.

[0122] Specifically, the method in step S2 includes:

[0123] S21. Define the cross product function: Let the starting point be O, then the vector from O to point A is denoted as... The vector pointing from point O to point B is denoted as ,definition and The cross product function, denoted as cross(), takes the coordinates of points O, A, and B as input parameters and returns the cross product value of a two-dimensional vector.

[0124] S22. Define a function to compute the convex hull of a point set, denoted as convex_hull(), with the point set points as input. First, sort the points in points in ascending order by their x-coordinates (or y-coordinates if x-coordinates are the same). Then, construct the convex hull in two steps: first, compute the lower convex hull by forward traversal, and then compute the upper convex hull by backward traversal. During construction, use the cross product (cross()) function to determine the direction of the points. If a new point causes the path to turn to the right or become collinear, continuously remove vertices from the stack to maintain convexity. Finally, remove duplicate endpoints from the upper and lower convex hulls and concatenate them to return the complete convex hull vertex sequence.

[0125] S23. Define a function for the minimum bounding rectangle, denoted as min_bounding_rectangle(), with the convex hull points as input parameters. Let the minimum area of ​​the rectangle be min_area, initialized to infinity.

[0126] Traverse each edge of the convex hull and obtain the orientation angle of the current edge, denoted as `angle`. Then, rotate the entire `hull_points` point set in the reverse direction by this angle, making the current edge parallel to the X-axis. The resulting new point set is denoted as `rotated_points`. In the rotated coordinate system, calculate the axis-aligned bounding box (i.e., the minimum and maximum x and y coordinates) of all points in `rotated_points`, thus obtaining the area of ​​the rectangle facing downwards. If this area is less than the recorded minimum value `min_area`, update the optimal solution: first, determine the four corner points of the bounding box at this point, then restore it to the original coordinate system through a forward rotation, saving it as the current best rectangle. Finally, after traversing all edges, return the vertex sequence of the rectangle with the smallest area as the minimum bounding box.

[0127] S24. Define a function to calculate the minimum bounding rectangle and the angles of adjacent sides, denoted as calc_mbr_and_angles(), with the input parameter being the polygon point set polygon_points.

[0128] See Figure 3 First, the convex hull of the polygon points is calculated using the convex_hull() function, denoted as hull. Then, the minimum bounding rectangle of hull is calculated using the min_bounding_rectangle() function, denoted as mbr. Next, with east as 0 degrees (i.e., the x-axis) and increasing counterclockwise, the angles of two adjacent sides of mbr are calculated, denoted as angle1 and angle2 respectively. Finally, mbr, angle1, and angle2 are returned as dictionaries. Figure 3 The effect of generating a minimum bounding rectangle for a polygon is shown, with the direction angles of the two adjacent sides added.

[0129] The core Python code for this process is as follows:

[0130] # Calculate the convex hull of a point set

[0131] hull = convex_hull(polygon_points)

[0132] # Calculate the minimum bounding rectangle

[0133] mbr = min_bounding_rectangle(hull)

[0134] # Calculate the angle between two adjacent sides (with due east as 0 degrees, increasing counterclockwise).

[0135] edge1 = (mbr[1][0] - mbr[0][0], mbr[1][1] - mbr[0][1]) # First neighboring edge vector

[0136] edge2 = (mbr[3][0] - mbr[0][0], mbr[3][1] - mbr[0][1]) # Second neighboring edge vector

[0137] # Convert angles to degrees

[0138] angle1 = math.degrees(math.atan2(edge1[1], edge1[0]))

[0139] angle2 = math.degrees(math.atan2(edge2[1], edge2[0]))

[0140] return {

[0141] 'mbr': mbr, # The four corner points of the minimum bounding rectangle

[0142] 'angle1': angle1, # Angle of the first adjacent edge (0 degrees for due east, counterclockwise)

[0143] 'angle2': angle2 # The angle of the second adjacent side (0 degrees for due east, counterclockwise)

[0144] }

[0145] S3. Add the rotation angle of the minimum rectangular bounding box to the GIS feature.

[0146] Specifically, the method in step S3 includes:

[0147] S31. Prepare the GIS vector data for which the rotation angle will be added, denoted as shp_file. Set the variable for the field to be added to output_field, with an initial value of 'rangle'.

[0148] S32. Import the Arcpy library and check if the `output_field` field name exists in the `shp_file`. If it does not exist, use the `arcpy.AddField_management()` function to add the `output_field` field to the `shp_file` with the type "FLOAT". If it exists, update the value of `output_field` to 'rotate_angle' before adding it to the `shp_file`.

[0149] S33. Use the arcpy.da.UpdateCursor() function to create an update cursor for the shp_file file based on the field list ["FID", "SHAPE@",output_field], denoted as cursor.

[0150] (1) Traverse the cursor, record the current record as row, and obtain the feature number fid (i.e. row[0]) and polygon object polygon (i.e. row[1]) based on row.

[0151] (2) Based on polygon, use traversal to obtain the set of vertices corresponding to the polygon, denoted as polygon_points.

[0152] S34. Use the calculate_mbr_and_angles() function to calculate the minimum bounding rectangle and adjacent side angles corresponding to the point set polygon_points, denoted as result.

[0153] Obtain the adjacent side angles of the minimum bounding rectangle according to the dictionary keywords "angle1" and "angle2", denoted as rangle1 = result['angle1'] and rangle2 = result['angle2'].

[0154] If abs(rangle1) < abs(rangle2), update rangle to rangle1; otherwise, update rangle to rangle2.

[0155] Update the attribute value for the current record row. Assign float(rangle) to row[2], and use cursor.updateRow(row) to update the record.

[0156] Refer to Figure 4 , after traversing all cursors, the operation of adding the rotation angle of the minimum rectangle bounding box to each feature of the shp_file file ends. Figure 4 Shows the effect of adding the rotation angle field "range" and assignment in the attribute table of the shp_file file.

[0157] S4. Create a texture mapping CGA rule in CityEngine.

[0158] Specifically, the method of step S4 includes:

[0159] S41. Create a CGA rule file topSkin.cga in CityEngine. In the file, use the attr keyword to declare the attribute Id as the texture index, with an initial value set to 0; at the same time, define the img() function to implement the function of dynamically returning the corresponding image file name according to Id.

[0160] S42. Use the attr keyword to declare the attribute rangle as the texture rotation angle (default value is 0); at the same time, declare the boolean attribute isRotate to control texture rotation (default value is false), and explicitly define its options as true and false through the @Enum annotation.

[0161] Define the calcRotateAngle() function, which makes a conditional judgment based on the isRotate attribute: if isRotate is true, the calculation result is 180 - rangle; if isRotate is false, the calculation result is -rangle.

[0162] S43. Use the `attr` keyword to define Boolean properties that control axial (X, Y, and Z) flipping: `xFlip`, `yFlip`, and `zFlip`. The default configuration is: flipping is enabled on the X-axis (true), and flipping is disabled on the Y and Z axes (false).

[0163] S44. Create the TopSkin top surface shape texture mapping rule.

[0164] S441. Adjust the model coordinate system. First, use the alignScopeToAxes(y) operation to align the current shape's range coordinate system to the global Y-axis. Then, use the alignScopeToGeometryBBox(xz) operation to align it to the XZ plane of the target geometry's bounding box. Next, execute the rotateScope(0, calcRotateAngle, 0) operation to rotate the coordinate system around the Y-axis by the angle calcRotateAngle. Finally, use the mirrorScope(xFlip, yFlip, zFlip) operation to mirror the model.

[0165] S442. Generate a rectangle based on the original shape. Use the primitiveQuad() operation to generate a quadrilateral that matches the bounding box of the original shape.

[0166] S443. Apply texture to the rectangle. The setupProjection(0,scope.xz,scope.sx,scope.sz) and projectUV(0) operations are used sequentially to set the texture layer and create the projected coordinate system for the rectangle. Finally, Skin is called to perform the texture application. Specifically, the Skin rule uses the fileExists() function to first check if the image img(Id) exists. If it does, the texture(img(Id)) operation is performed; otherwise, the NIL operation is used to delete the current shape.

[0167] The core code for the above process is as follows:

[0168] # Filename: topSkin.cga

[0169] @Hidden

[0170] attr Id=0 # Image ID

[0171] img(Id)="assets / png / "+Id+".png"

[0172] @Hidden

[0173] attr rangle=0 # Texture rotation angle

[0174] @Description("Whether to rotate the texture")

[0175] @Enum(true, false)

[0176] attr isRotate=false

[0177] # Define rotation angle function

[0178] calcRotateAngle=case isRotate: 180-(rangle) else: -(rangle)

[0179] @Description("X-axis flip")

[0180] attr xFlip = true

[0181] @Description("Y-axis flip")

[0182] attr yFlip = false

[0183] @Description("Z-axis flip")

[0184] attr zFlip = false

[0185] # Texture mapping rules

[0186] TopSkin -->

[0187] alignScopeToAxes(y)

[0188] alignScopeToGeometryBBox(xz)

[0189] rotateScope(0,calcRotateAngle,0)

[0190] mirrorScope(xFlip, yFlip, zFlip)

[0191] primitiveQuad()

[0192] setupProjection(0,scope.xz,scope.sx,scope.sz)

[0193] projectUV(0)

[0194] Skin # Execution Texture

[0195] Skin --> # Apply texture to an image based on the specified ID

[0196] case fileExists(img(Id)):

[0197] texture(img(Id))

[0198] else: NIL

[0199] S5. Apply texture mapping and correct the irregular top surface of the 3D model.

[0200] Specifically, the method in step S5 includes:

[0201] S51. Import the GIS vector data `shp_file` into CityEngine to generate shape layers. Verify that the layers contain `Id` and `rangle` attributes, and generate 3D models for each shape (i.e., the features in `shp_file`). Then, select the top face shape of all models and apply the `TopSkin` rule from the `topSkin.cga` file to generate a model with preliminary texture (denoted as the rough model), the effect of which is shown below. Figure 5 As shown.

[0202] S52. Manually correct textures that do not match the shape. Although the rough model generated through the preceding steps can achieve correct texture mapping for most target shapes, there are still a few targets where the texture orientation does not match the geometry (i.e., the texture orientation is abnormal). To address this, manual adjustment is performed using the inspector. Specifically, adjust the Boolean property controlling texture rotation, `isRotate`, to `true`, rotating the texture by 180 degrees from its initial angle to ensure precise alignment with the current shape. The corrected texture mapping effect is as follows. Figure 6 As shown, the final effect is as follows Figure 7 As shown.

[0203] Through the above steps, this invention achieves fast and accurate texturing of the top surface shape of a 3D model (regardless of whether it is regular or not). This method significantly improves texturing efficiency while ensuring texturing quality.

[0204] In large-scale 3D modeling of geographic information systems, this application enables batch texturing of the top surface of 3D models. Based on GIS vector and raster data, the generated texture images strictly adhere to the shape of the top surface of each model (regardless of whether it is regular or not), thus ensuring the accuracy of the texture files. Furthermore, by calculating and adjusting the texture rotation angle, texturing accuracy is further guaranteed. As an automated batch processing solution, compared to traditional manual texturing, this application significantly saves labor costs and improves operational efficiency.

[0205] The above description is merely an embodiment of this application and does not limit the patent scope of this application. Any equivalent structural or procedural transformations made using the content of this application's specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

Claims

1. A batch texture mapping method for irregular top surfaces of 3D models, characterized in that, include: S1. Generate matching texture images based on GIS features; S2. Calculate the rotation angle of the smallest rectangular bounding box of the GIS feature; S3. Add the rotation angle of the minimum rectangular bounding box to the GIS element; S4. Create a texture mapping CGA rule in CityEngine; S5. Apply texture mapping and correct the irregular top surface of the 3D model.

2. The method according to claim 1, characterized in that, The method of step S1 includes: S11. Prepare GIS vector data shp_file and raster data tif_file that are fully corresponding to the coordinate system, configure the texture image output folder output_dir, the field name id_field used as feature identifier in the shp_file, and the invalid raster value no_data_value in the tif_file, with the default value of no_data_value being 0; S12. Import the Arcpy and PIL.Image libraries, set arcpy.env.overwriteOutput = True, and check if output_dir, shp_file, tif_file and id_field exist in turn. If they do not exist, throw an exception. S13. Read the id_field and feature geometry object field "SHAPE@" in the shp_file using the arcpy.da.SearchCursor() function, establish a cursor and convert it into a feature list, traverse the feature list, and obtain the feature_id and feature_shape of each feature; S14. Call the arcpy.management.MakeRasterLayer() function to register the tif_file as a temporary in-memory raster map raster_layer, and use the arcpy.management.Clip() function to clip the raster_layer according to the feature_shape to obtain a temporary clipped raster file temp_raster; S15. Use the arcpy.RasterToNumPyArray() function to convert the temp_raster into a NumPy array raster_array. Based on the dimension of the raster_array, perform classification processing to generate an RGBA format PNG texture image with an alpha channel. If the image is invalid, output a prompt message. After traversing all elements, delete the temp_raster.

3. The method according to claim 2, characterized in that, The specific method for classifying based on the raster_array dimension in step S15 is as follows: If the raster_array has a dimension of 3 and the number of bands is not less than 3, it is a multi-band image. The first three bands are extracted as red, green and blue channels, respectively. Each channel is normalized and converted to uint8 type. They are then merged into rgb_array using np.dstack. An alpha channel is set for rgb_array according to no_data_value. The rgb_array and alpha_channel are then merged into rgba_array again using np.dstack. The image is generated in 'RGBA' mode using the Image.fromarray() function and saved as a png. If the raster_array has a dimension of 2, it is a single-band image. The raster_array is normalized, and an alpha channel is created based on no_data_value using the np.where() function. The normalized data and the alpha channel are merged into an rgba_array. The image is generated in 'RGBA' mode using the Image.fromarray() function and saved as a png.

4. The method according to claim 1, characterized in that, The method of step S2 includes: S21. Define the cross product function cross(), which takes the coordinates of the starting point O, point A, and point B as input and returns the cross product value of the two-dimensional vectors pointing from O to A and from O to B. S22. Define a convex hull function convex_hull(), take a point set points as input, sort the point set in ascending order of x-coordinate, and sort in ascending order of y-coordinate if the x-coordinates are the same, determine the direction of the points through the cross product function, construct the lower convex hull and the upper convex hull respectively, remove duplicate endpoints and concatenate to obtain the convex hull vertex sequence and return it. S23. Define the minimum bounding box function min_bounding_rectangle(), take the convex hull point set hull_points as input, traverse each edge of the convex hull, calculate the direction angle of the current edge and rotate the hull_points in the reverse direction until the current edge is parallel to the X-axis, calculate the area of ​​the axis-aligned bounding box of the rotated point set, filter out the bounding box with the smallest area, restore its corner points to the original coordinate system, and return the vertex sequence of the minimum bounding box; S24. Define a calculation function calc_mbr_and_angles(), take the polygon point set polygon_points as input, call the convex hull function and the minimum bounding box function in sequence to obtain the minimum bounding rectangle mbr of the feature, calculate the direction angles angle1 and angle2 of two adjacent sides of the mbr with due east as 0 degrees and counterclockwise angles increasing as the rule, and return mbr, angle1 and angle2 in dictionary form.

5. The method according to claim 4, characterized in that, The specific rules for determining the direction of a point using the cross product function in step S22 are as follows: If the cross product is less than 0 after a new point is added, it indicates that the path turns to the right, and the vertex is removed from the stack; if the cross product is equal to 0, it indicates that the points are collinear, and the vertex is removed from the stack. If the cross product is greater than 0, it indicates that the path turns to the left. The new point is retained and added to the stack.

6. The method according to claim 4, characterized in that, The specific steps for calculating the directional angles of two adjacent sides of the mbr in step S24 are as follows: Extract the coordinates of the adjacent vertices of the MBR, and calculate the vectors edge1 and edge2 of the two adjacent edges; calculate the radian values ​​of edge1 and edge2 respectively using the math.atan2() function, and then convert the radian values ​​to angle values ​​using the math.degrees() function to obtain angle1 and angle2.

7. The method according to claim 1, characterized in that, The method of step S3 includes: S31. Prepare the GIS vector data shp_file to which the rotation angle will be added, and set the initial name of the field to be added to 'rangle'; S32. Import the Arcpy library, check if the 'rangle' field exists in the shp_file. If it does not exist, add a 'rangle' field of type FLOAT using the arcpy.AddField_management() function. If it exists, update the field name to 'rotate_angle' and add a field of the same type. Denote the newly added field as output_field. S33. Read the "FID", "SHAPE@" and output_field fields of the shp_file through the arcpy.da.UpdateCursor() function, create an update cursor and traverse it to obtain the feature number fid and the polygon object polygon, and extract the vertex set polygon_points of the polygon; S34. Call the calc_mbr_and_angles() function to obtain angle1 and angle2, select the angle with the smaller absolute value as the feature rotation angle rangle, assign float(rangle) to the output_field field, and update the feature attribute record through cursor.updateRow(row). After traversing all cursors, the rotation angle is added.

8. The method according to claim 1, characterized in that, The method of step S4 includes: S41. Create a CGA rule file topSkin.cga, use the attr keyword to declare the attribute Id as the texture index, initialize the value to 0, and define the img() function to dynamically return the file path and name of the corresponding texture image based on the Id; S42. Use the `attr` keyword to declare the attribute `rangle` as the texture rotation angle, with an initial value of 0; declare the boolean attribute `isRotate` to control the secondary rotation of the texture, with an initial value of false, and explicitly define its optional values ​​as true and false using the `@Enum` annotation; define the function `calcRotateAngle()`, which returns 180 - `rangle` if `isRotate` is true, and returns - `rangle` if `isRotate` is false. S43. Use the attr keyword to declare the boolean properties xFlip, yFlip, and zFlip that control the mirror flipping of the X, Y, and Z axes. The default configuration is xFlip=true, yFlip=false, and zFlip=false. S44. Create the TopSkin core rule for top surface texture mapping, and sequentially perform coordinate system alignment, rotation, and mirroring to generate a quadrilateral consistent with the original shape bounding box of the feature. Set the texture projection layer and coordinate system, and call the Skin rule to perform the mapping operation. The Skin rule first checks whether the texture image corresponding to img(Id) exists through the fileExists() function. If it exists, it performs the texture(img(Id)) mapping operation; otherwise, it performs the NIL operation to delete the current shape.

9. The method according to claim 8, characterized in that, In step S44, The specific CGA operation commands for performing coordinate system alignment, rotation, and mirroring are as follows: First, alignScopeToAxes(y) is executed to align the current shape's coordinate system to the global Y-axis. Then, alignScopeToGeometryBBox(xz) is executed to align the coordinate system to the XZ plane of the target geometry's bounding box. Next, rotateScope(0, calcRotateAngle, 0) is executed to rotate the coordinate system around the Y-axis by the angle calcRotateAngle. Finally, mirrorScope(xFlip, yFlip, zFlip) is executed to complete the coordinate system mirroring. The specific CGA operation commands for generating quadrilaterals and setting texture projection are as follows: The primitiveQuad() operation generates a quadrilateral that matches the bounding box for the original shape; the setupProjection(0, scope.xz, scope.sx, scope.sz) and projectUV(0) operations are executed in sequence to complete the texture layer setup and projection coordinate system creation.

10. The method according to claim 1, characterized in that, The method of step S5 includes: S51. Import the GIS vector data shp_file with Id and rangle attributes into CityEngine, generate a shape layer and verify the integrity of the attributes, generate a 3D model based on the shape layer, select the top surface shape of all models, apply the TopSkin rule, and generate a rough model with preliminary texture. S52. For models in the coarse model where the texture direction does not match the geometry, adjust the isRotate property to true in the attribute inspector, so that the texture is rotated 180 degrees on the basis of the initial rotation angle, to complete the precise correction of the texture and realize the batch texture mapping of the irregular top surface of the 3D model.

Images