Light intensity distribution control method for volumetric bioprinting based on projection algorithm

Abstract

The present invention proposes a volume bioprinting light intensity distribution control method based on a projection algorithm, comprising: slicing a three-dimensional model equidistantly along the Z axis to obtain M transverse slices; performing a 360-degree Laden transform on each transverse slice; stacking filtered projection data of the same angle of the M transverse slices in sequence along the direction of the Z axis to form a projection slice of the same angle; taking any M1 of the M transverse slices as transverse slice samples and executing a volume bioprinting evaluation algorithm to obtain an optimal window function; filtering the projection slices using the optimal window function and based on a filtered back-projection algorithm; performing a printing light intensity adjustment operation on the filtered projection slices to obtain pre-projection slices; performing volume bioprinting using the pre-projection slices, and projecting the pre-projection slices in sequence according to angle order into a printing bottle containing the specific bio-ink.

Description

Technical Field

[0001] The present invention belongs to the technical field of 3D printing, and in particular relates to a volumetric bioprinting light intensity distribution control method based on a projection algorithm. Background Art

[0002] At present, the conventional technology of volumetric bioprinting is to use layer-by-layer printing technology. The so-called layer-by-layer printing is to slice the three-dimensional model perpendicular to the Z axis. During printing, these XY axis two-dimensional slices are printed vertically in the Z axis direction in the order of their positions. Finally, a three-dimensional object is stacked layer by layer from "two-dimensional" to "three-dimensional".

[0003] There are currently two main implementation schemes for this layer-by-layer printing technology, namely extrusion and light-curing. The extrusion method uses air pressure or mechanically driven nozzles to controllably extrude the bio-ink. The (bio) ink is extruded from the nozzle and deposited on the forming platform to form a two-dimensional structure. As the nozzle or the forming platform moves in the z direction, the two-dimensional structure accumulates layer by layer to form a three-dimensional structure. Light-curing printing uses a digital light projector to cure the entire surface of the bio-ink, and through the up and down movement of the forming platform, it is cured layer by layer to obtain a three-dimensional structure.

[0004] However, with this layer-by-layer printing technology, the higher the printing accuracy required, the faster the number of two-dimensional slices cut out, resulting in a long printing time, generally from tens of minutes to several hours. And too long a wait will cause a large number of living cells in the biological ink to die. Secondly, because each slice is formed separately in sequence during layer-by-layer printing, the cells and vascular networks in the printed object cannot be organically integrated, and there will be a height difference between the layers on the surface, which will cause the printed object to form a layered rough texture on the surface, which is inconsistent with the relatively smooth surface of normal organs. Finally, in order to reduce the thickness of the printed slices, the light-curing layer-by-layer printing technology uses weakly penetrating ultraviolet light to irradiate the biological ink for curing, which is extremely unfriendly to living cells and will further reduce cell survival rates. Summary of the invention

[0005] The purpose of the present invention is to propose a volume bioprinting light intensity distribution control method, device, equipment and storage medium based on projection algorithm, which can effectively solve the above-mentioned technical problems existing in the prior art.

[0006] In order to achieve the above object, an embodiment of the present invention provides a method for controlling light intensity distribution of volumetric bioprinting based on a projection algorithm, the method comprising the steps of:

[0007] S1, slicing the three-dimensional model horizontally along the Z axis at equal intervals to obtain M horizontal slices; wherein each horizontal slice obtained is a horizontal slice binary image, M≥1;

[0008] S2, performing a 360-degree Laden transformation on each of the transverse slices, so that each transverse slice obtains W lines of projection data; wherein each line of projection data is projection data of a Laden transformation at an angle of each of the transverse slices, and among the W lines of projection data obtained by performing the Laden transformation on the same transverse slice, the wth line of projection data is projection data of a Laden transformation at an angle (w-1)*θ of the transverse slice, w=1, 2...W, 0°<θ≤1°, W*θ=360°;

[0009] S3, stacking the projection data of the same angle of the M transverse slices in sequence along the direction of the Z axis to form a projection slice of the same angle, thereby obtaining projection slices of the side of the three-dimensional model at different angles of 360 degrees;

[0010] S4, taking any M1 of the M transverse slices as transverse slice samples and executing a volumetric bioprinting evaluation algorithm to obtain an optimal window function; wherein the optimal window function has a corresponding critical curing brightness value;

[0011] S5, using the optimal window function and based on the filtered back projection algorithm, filtering each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model, so as to obtain filtered projection slices;

[0012] S6, performing a printing light intensity adjustment operation on the filtered projection slice, by simultaneously adjusting the printer light intensity and adjusting the pixel value of the filtered projection slice so that the specific bio-ink is cured when the pixel value of the projection slice is greater than or equal to the critical curing brightness value, thereby obtaining a pre-projection slice;

[0013] Wherein, the step S6 comprises:

[0014] S61, adjusting the light intensity of the projector of the printer so that the light intensity projected by the projector when the pixel value is close to the critical curing brightness value is equal to the critical curing light intensity value of the specific biological ink;

[0015] S62, adjusting the pixel values ​​of all pixel points of the filtered projection slice by the following pixel value adjustment algorithm, so that the pixel value of the pixel point corresponding to the critical solidification brightness value is adjusted to the adjusted pixel value, and using the projection slice after the pixel value adjustment as the pre-projection slice:

[0016]

[0017] Among them, h is the adjusted pixel value, b is the pixel value adjustment parameter, h o is the original pixel value, <> represents rounding;

[0018] S7. Apply the pre-projected slices to perform volumetric bio-printing, project the pre-projected slices sequentially into a printing bottle containing the specific bio-ink in angular order, and rotate the printing bottle at a constant speed so that the angle of the projected slices corresponds to the angle of the printing bottle to perform volumetric printing of the three-dimensional model; wherein the bio-inks of all transverse sections in the printing bottle are simultaneously photocured into corresponding transverse slices of the three-dimensional object.

[0019] As an improvement of the above solution, step S7 includes:

[0020] S71, start the stepper motor, control the stepper motor to rotate at a constant speed of K degrees per second, so as to drive the printing bottle arranged on the stepper motor to rotate synchronously; wherein the printing bottle contains the specific biological ink, 0 <K;

[0021] S72, start the projection device to project the pre-projected slice onto the side of the printing bottle, so that the position of the pre-projected slice on the projection screen satisfies:

[0022] X=[(x 0 -x)÷2]

[0023] Y=[(y 0 -y)÷2]

[0024] Wherein, (X, Y) is the coordinate of the upper left corner of the pre-projected slice on the projection screen, and the coordinate is calculated by taking the upper left corner of the projection screen as the origin (0, 0), the horizontal rightward direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; 0 ,y 0 is the width and height of the projection screen; x, y is the width and height of the pre-projected slice;

[0025] S73, controlling the pre-projected slices to be projected onto the side of the printing bottle in angular order, and the projection speed is the same as that of the stepping motor, so that the angle between the pre-projected slices and the printing bottle always remains the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

[0026] As an improvement of the above solution, a print size adjustment step is further included between step S5 and step S6, and the print size adjustment step includes:

[0027] Performing a border removal operation on each of the filtered projection slices to obtain a valid projection slice; the border removal operation includes removing all black columns on the left and right sides, removing all black rows on the top and bottom sides, and retaining only the valid area in the middle;

[0028] According to the magnification or reduction factor, multiplying the magnification or reduction factor by the side length of the effective projection slice to obtain the scaled projection slice size;

[0029] Using an image scaling algorithm to scale the effective projection slice according to the magnification or reduction multiple, when scaling up, if the size of the scaled projection slice exceeds the size of the projection screen, then calculating the portion of the middle effective area that has the same size as the projection screen, and using the calculated projection slice as the scaled projection slice;

[0030] The filtered projection slice in step S6 is the scaled projection slice.

[0031] As an improvement of the above solution, the image scaling algorithm is a bi-triple interpolation algorithm.

[0032] As an improvement of the above scheme, step S7 also includes a printing time control step, in which the timing starts from when the projection device is started to start projection, and printing is stopped when the printing time parameter is reached; wherein, when stopping printing, the projection is first turned off, and then the rotation of the stepper motor is stopped.

[0033] As an improvement of the above solution, step S4 specifically includes:

[0034] S41, taking M1 transverse slices obtained by equidistant transverse slicing of the three-dimensional model along the Z axis as transverse slice samples; wherein each transverse slice sample is a transverse slice binary image, and M1≥1;

[0035] S42, performing a 360-degree Laden transformation on each of the transverse slice samples, so that each of the transverse slice samples obtains W lines of projection data; wherein each line of projection data is projection data of an angle of each of the transverse slice samples obtained by Laden transformation, and among the W lines of projection data obtained by performing Laden transformation on the same transverse slice sample, the wth line of projection data is projection data of an angle of (w-1)*θ of the transverse slice sample obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°;

[0036] S43, performing filtering back-projection transformation on the W lines of projection data obtained from each transverse slice sample using a to-be-determined window function, thereby obtaining M1 back-projection reconstruction images with the same resolution as each transverse slice sample; wherein the back-projection reconstruction image is the cumulative distribution of light intensity of the transverse slice at the corresponding height of the printed object in the volume bioprinting projection;

[0037] S44, normalizing the pixel value of each back-projection reconstruction image to an integer between 0 and 255 to obtain M1 normalized reconstruction images;

[0038] S45, setting a variable parameter critical curing brightness value, respectively setting the critical curing brightness value to different integers of [1, 254], and setting the uncured pixels in each normalized reconstructed image to 0 and the cured pixels to 1 according to the critical curing brightness value, thereby correspondingly obtaining M1*254 simulated printed transverse slices; wherein, the pixels in the normalized reconstructed image that are less than the critical curing brightness value are uncured pixels, and the pixels that are greater than or equal to the critical curing brightness value are cured pixels;

[0039] S46, calculating the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height by the following formula (1), to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height:

[0040]

[0041] Where P is the similarity evaluation value, X and Y are the horizontal resolution and vertical resolution of the horizontal slice sample and the simulated printed horizontal slice; f o (x, y) is the value of the pixel at the position (x, y) in the horizontal slice sample with the upper left corner as the coordinate (1, 1) and the positive direction as downward and left. f(x, y) is the value of the pixel at the position (x, y) in the horizontal slice of the simulated print with the upper left corner as the coordinate (1, 1) and the positive direction as downward and left.

[0042] S47, respectively calculating the average values ​​of M1 similarity evaluation values ​​corresponding to M1 simulated printed transverse slices at the same critical curing brightness value, and taking the average value closest to 0 as the printing evaluation index of the to-be-determined window function;

[0043] S48, using different window functions to perform filtered back projection transformation on the W lines of projection data obtained for each of the transverse slice samples in step S43, and comparing the printing evaluation indicators of the different window functions obtained after steps S44 to S47, and taking the window function corresponding to the printing evaluation indicator closest to 0 among the printing evaluation indicators of different window functions as the optimal window function.

[0044] As an improvement of the above solution, in step S46, the following formula (2) is used instead of formula (1) to calculate the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height:

[0045]

[0046] Wherein, in said S47, the average value closest to 1 is taken as the printing evaluation index of the undetermined window function; in said S48, the window function corresponding to the printing evaluation index closest to 1 among the printing evaluation indexes with different parameters is taken as the optimal window function.

[0047] As an improvement of the above solution, step S5 specifically includes:

[0048] S51, performing a fast Fourier transform on each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model; wherein the width N used in the fast Fourier transform is greater than or equal to the number of pixels of each line of projection data of each transverse slice;

[0049] S52, using the optimal window function to perform a windowing operation on the ramp filter, and using the ramp filter after the windowing operation to filter each line of projection data of each projection slice after fast Fourier transform;

[0050] S53, performing inverse fast Fourier transform on each line of projection data of each projection slice after filtering, so as to obtain a filtered projection slice.

[0051] As an improvement of the above solution, step S52 specifically includes:

[0052] S521, discretizing the ramp filter in the frequency domain for sampling; the number of sampling points is an even number greater than or equal to N and closest to N;

[0053] S522, performing point-to-point multiplication of the discretized ramp filter and the optimal window function with the same number of discretized sampling points to perform a windowing operation;

[0054] S523, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head;

[0055] S524, perform point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after the fast Fourier transform of each projection slice to perform a filtering operation, and discard the extra points of the ramp filter after the windowing and translation operation compared with the points of each line of projection data after the fast Fourier transform of each projection slice.

[0056] As an improvement of the above solution, in step S2, Laden transform is performed on each angle of 360 degrees on the side of each transverse slice in turn, so as to obtain Laden transformed projection data of the corresponding angle of each transverse slice.

[0057] Compared with the prior art, the volume bioprinting light intensity distribution control method based on projection algorithm provided by the embodiment of the present invention has the following specific technical effects:

[0058] 1. Faster printing speed (10 to 120 seconds)

[0059] ① Printing with a limited number of filtered projection slices, which greatly reduces the number of slices compared to conventional layer-by-layer printing.

[0060] 2. Higher cell survival rate (more than 95%)

[0061] ①. Due to the fast printing speed, cell damage caused by long-term printing is avoided.

[0062] ②. There is no direct contact with the printed object during the printing process, avoiding cell damage and contamination.

[0063] ③. The printing process is kept at room temperature and uses visible light to print, avoiding damage to cells caused by temperature, laser or ultraviolet light.

[0064] 3. The printing target is formed in one piece, solving the problem that cells and vascular networks cannot be organically integrated.

[0065] ① Each filtered projection slice is obtained by filtering the Laden transform data of the same angle of all transverse slices and stacking them together. Therefore, when printing, all transverse slices are photocured from the side at the same time, avoiding the problem of the inability of cells and vascular networks to be organically integrated due to layer-by-layer printing. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] In order to more clearly illustrate the technical solution of the present invention, the drawings required for use in the implementation mode will be briefly introduced below. Obviously, the drawings described below are only some implementation modes of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0067] Figure 1 A flow chart of a method for controlling light intensity distribution in volumetric bioprinting based on a projection algorithm provided in an embodiment of the present invention.

[0068] Figure 2 A schematic structural diagram of a three-dimensional model (bone screw 3D model) of a volumetric bioprinting light intensity distribution control method based on a projection algorithm provided for implementing an embodiment of the present invention.

[0069] Figure 3 A schematic diagram of the Radon transform performed on a transverse slice sample by a volume bioprinting light intensity distribution control method based on a projection algorithm provided in an embodiment of the present invention.

[0070] Figure 4 Demonstrate the use of a volume bioprinting light intensity distribution control method based on a projection algorithm provided by an embodiment of the present invention Figure 2 The three-dimensional model shown is a schematic diagram of the entire printing process.

[0071] Figure 5a The process of generating (transforming) the obtained projection slices using a volume bioprinting control method provided by an embodiment of the present invention is demonstrated.

[0072] Figure 5b A schematic diagram of volumetric printing of a three-dimensional model using a projection slice obtained by using a volumetric bio-printing light intensity distribution control method based on a projection algorithm provided in an embodiment of the present invention.

[0073] Figure 6 yes Figure 2 A schematic diagram of a 0 degree unfiltered projection slice of the bone screw 3D model is shown.

[0074] Figure 7 This is a line graph showing the average evaluation value of the Kaiser window function with different β values ​​changing with the critical curing light intensity value.

[0075] Figure 8 It is a line graph of the optimal evaluation value of the Kaiser window function with different β values.

[0076] Fig. 9 It is a line graph showing the average evaluation value of the Kaiser window function with a β of 5 and the change in the critical curing light intensity value.

[0077] Fig.10 yes Figure 2 The bone screw 3D model shown is a schematic diagram of a 0-degree filtered projection slice after being filtered using a ramp filter windowed with a Kaser window function of β=5.

[0078] Fig.11 yes Figure 2 The bone screw 3D model shown is a schematic diagram of a scaled projection slice after filtering at 0 degrees and removing invalid parts using a ramp filter with a Kaser window function of β=5. DETAILED DESCRIPTION

[0079] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0080] Figure 1A flow chart of a method for controlling light intensity distribution of volumetric bioprinting based on a projection algorithm provided in an embodiment of the present application, the method comprising steps S1 to S7:

[0081] S1, slicing the three-dimensional model horizontally along the Z axis at equal intervals to obtain M horizontal slices; wherein each horizontal slice obtained is a horizontal slice binary image, M≥1;

[0082] S2, performing a 360-degree Laden transformation on each of the transverse slices, so that each transverse slice obtains W lines of projection data; wherein each line of projection data is projection data of a Laden transformation at an angle of each of the transverse slices, and among the W lines of projection data obtained by performing the Laden transformation on the same transverse slice, the wth line of projection data is projection data of a Laden transformation at an angle (w-1)*θ of the transverse slice, w=1, 2...W, 0°<θ≤1°, W*θ=360°;

[0083] S3, stacking the projection data of the same angle of the M transverse slices in sequence along the direction of the Z axis to form a projection slice of the same angle, thereby obtaining projection slices of the side of the three-dimensional model at different angles of 360 degrees;

[0084] S4, taking any M1 of the M transverse slices as transverse slice samples and executing a volumetric bioprinting evaluation algorithm to obtain an optimal window function; wherein the optimal window function has a corresponding critical curing brightness value;

[0085] S5, using the optimal window function and based on the filtered back projection algorithm, filtering each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model, so as to obtain filtered projection slices;

[0086] S6, performing a printing light intensity adjustment operation on the filtered projection slice, by simultaneously adjusting the printer light intensity and adjusting the pixel value of the filtered projection slice so that the specific bio-ink is cured when the pixel value of the projection slice is greater than or equal to the critical curing brightness value, thereby obtaining a pre-projection slice;

[0087] S7. Apply the pre-projected slices to perform volumetric bio-printing, project the pre-projected slices sequentially into a printing bottle containing the specific bio-ink in angular order, and rotate the printing bottle at a constant speed so that the angle of the projected slices corresponds to the angle of the printing bottle to perform volumetric printing of the three-dimensional model; wherein the bio-inks of all transverse sections in the printing bottle are simultaneously photocured into corresponding transverse slices of the three-dimensional object.

[0088] Wherein, the step S6 comprises:

[0089] S61, adjusting the light intensity of the projector of the printer so that the light intensity projected by the projector when the pixel value is close to the critical curing brightness value is equal to the critical curing light intensity value of the specific biological ink;

[0090] S62, adjusting the pixel values ​​of all pixel points of the filtered projection slice by the following pixel value adjustment algorithm, so that the pixel value of the pixel point corresponding to the critical solidification brightness value is adjusted to the adjusted pixel value, and using the projection slice after the pixel value adjustment as the pre-projection slice:

[0091]

[0092] Among them, h is the adjusted pixel value, b is the pixel value adjustment parameter, h o is the original pixel value, and <> represents rounding.

[0093] Below, the various steps of a volume bioprinting light intensity distribution control method based on a projection algorithm provided in an embodiment of the present application will be described in detail in conjunction with the relevant drawings.

[0094] First, in the S1, the three-dimensional coordinates of a large number of triangles can be stored, and then the triangles that meet the conditions can be calculated, the coordinates of the triangle vertices can be calculated, and then the points can be connected to obtain a contour map, the contour map can be pixelated, and all the plane triangles can be moved up a fixed distance, and the same operation can be performed until the top layer of the three-dimensional model is reached or the top layer is exceeded for the first time in the plane parallel to the X and Y axes.

[0095] The smaller the fixed distance parallel to the X and Y axis planes is moved upward through step S1, the better. Because the number of horizontal slice binary images determines the accuracy of the three-dimensional object finally printed. The number of horizontal slice binary images is obtained according to the total height of the 3D model / fixed moving distance. In this embodiment, the number of horizontal slice binary images M≥1. It can be understood that equidistant horizontal slicing of the three-dimensional model along the Z axis is a discretization of it, which is bound to lead to a sufficient number of discrete faces. It is preferred that the three-dimensional model is more similar to the original three-dimensional model when the number of equidistant horizontal slices along the Z axis is sufficient. The printed object has a higher degree of restoration.

[0096] Specifically, the following is combined Figure 2 and Figure 4 , how to obtain the binary image of each transverse slice of the three-dimensional model will be described in detail. Figure 2 The bone screw 3D model shown is used as an example of the three-dimensional model of this embodiment. In S1, the three-dimensional model is composed of a large number of triangles, each triangle has a corresponding vertex three-dimensional coordinate, and the expression parallel to the X and Y axis plane is: Z = z 0 , where Z represents the coordinate of the Z axis, z0 is a constant, and z 0 Initialized to 0; the step S1 specifically includes steps S11 to S15:

[0097] S11. Compare the Z-axis coordinates of the three vertices of each triangle with the z 0 Compare and select the ones that satisfy the condition that the Z coordinate is greater than z 0 The sum of the vertices is less than z 0 The vertex of the triangle or two vertices have a Z coordinate equal to z 0 triangle.

[0098] S12, for the existence of Z coordinate greater than z 0 The sum of the vertices is less than z 0 The vertex of the triangle exists when the Z coordinate is equal to z 0 When there are vertices, connect the remaining two vertices and calculate the intersection of the connected line segment and the plane parallel to the X and Y axes.

[0099] S13, for the case where there is no Z coordinate equal to z 0 For the vertices of the triangle, calculate two Z coordinates greater than z 0 The vertices and the Z coordinates are less than z 0 The intersection of the line segment connecting the vertices and the plane parallel to the X and Y axes, or the intersection of two Z coordinates less than Z 0 The vertices and the Z coordinates are greater than z 0 The intersection of the line segment connecting the vertices and the plane parallel to the X and Y axes;

[0100] In the above steps S22 and S23, each intersection point is calculated using the vector method, and the formula is as follows:

[0101]

[0102]

[0103] Among them, the Z coordinate is greater than z 0 Point P 1 The coordinates of (x 1 ,y 1 , z 1 ), another Z coordinate is less than z 0 Point P 2 The coordinates of (x 2 ,y 2 , z 2 ), the intersection point is P, the plane expression parallel to the X and Y axes is ax+by+cz+d=0, a, b, c, d are constants, x, y, z are the variables of the X axis, Y axis, and Z axis, and the origin is O, is the vector from O to P, From O to P 1 The vector of P 1 To P 2 A vector of

[0104] Since the plane parallel to the X and Y axes is parallel to the X and Y axes, a = 0, b = 0, d = -z 0 , the formula is transformed into:

[0105]

[0106] in, Indicates P 1 The vector to P;

[0107] Thus, the coordinates of each intersection point P are obtained as Therefore, each triangle gets an intersection point, a Z coordinate equal to z 0 The endpoints or two intersection points or two Z coordinates equal to z 0 The endpoints of a triangle are taken as a group of two points.

[0108] S14, connect the two points of all point groups to obtain the Z-axis coordinate of the three-dimensional model as z 0 , and the contour image of the cross section parallel to the X and Y axis planes, the contour image is pixelated, the pixel points inside the contour image are white, and the pixel points outside the contour image are black, so as to obtain the horizontal slice binary image of the current layer. Among them, Figure 4 The transverse slice binary images (original binary images) of the 50th transverse slice and the 200th transverse slice of the three-dimensional model are respectively shown.

[0109] It can be understood that in this step, the contour map is pixelated into each pixel of the contour map, the center point coordinates are taken as the pixel coordinates, and the ray method is used to determine whether it is inside the contour. The pixel points inside the contour are white, and the pixel points outside the contour are black. The black represents that this pixel point is outside the three-dimensional model and is empty, indicating that the pixel point position does not need to be photocured into an object; the white represents that this pixel point is inside the three-dimensional model and is a solid, indicating that the pixel point position needs to be photocured into an object. Among them, the specific operation of pixelating the contour map is: the contour map is covered with square grids of the same size, each grid is a pixel, and each grid is set to a color, and the side length of the grid is a self-defined length.

[0110] S15, determine the current z 0 Is it less than H? If so, set z 0 +h value replaces the current z 0 And return to step S11; if not, end; wherein H is the total height of the three-dimensional model, h<H. As a preferred solution, h≤0.01mm.

[0111] It can be understood that in step S15, h is the fixed distance to move the plane parallel to the X and Y axes upward. The total height of the three-dimensional model is calculated by rounding. Specifically, the total height of the three-dimensional model / the rounded fixed distance moved upward = the number of layers + 1, that is, [H / h] = M+1. When the total height of the three-dimensional model is the same, the smaller the fixed distance moved upward, the more layers there are. The highest point of the three-dimensional model is located in the horizontal slice of the top layer.

[0112] It can be understood that in step S15, h is the fixed distance to move the plane parallel to the X and Y axes upward. The total height of the three-dimensional model is calculated by rounding. Specifically, the total height of the three-dimensional model / the rounded fixed distance moved upward = the number of layers + 1, that is, [H / h] = M+1. When the total height of the three-dimensional model is the same, the smaller the fixed distance moved upward, the more layers there are. The highest point of the three-dimensional model is located in the horizontal slice of the top layer.

[0113] It can be understood that after obtaining the binary images of the transverse slices of all transverse slices of the three-dimensional model through step S1, the embodiment of the present invention does not directly use the transverse slices as projection slices for 3D volume bioprinting, but performs a 360-degree Radon transform on all transverse slices of the three-dimensional model to obtain the Radon transformed projection data of different angles of all transverse slices, and stacks the projection data of the same angle of all transverse slices in sequence along the direction of the Z axis to form projection slices of the same angle, thereby obtaining projection slices of the side of the three-dimensional model at different angles of 360 degrees, and then uses the obtained projection slices to perform 3D volume bioprinting in angular order.

[0114] Specifically, in step S2, Radon transformation is performed on each angle of the 360-degree side of each transverse slice in sequence, so as to obtain Radon transformed projection data of the corresponding angle of each transverse slice. For example, taking θ=1° and N=360 as an example, Radon transformation is performed on the 0° angle of the 360-degree side of each transverse slice to obtain the first row of projection data, and then Radon transformation is performed on the 1° angle of the 360-degree side of each transverse slice to obtain the second row of projection data, and then Radon transformation is performed on the 2° angle of the 360-degree side of each transverse slice to obtain the third row of projection data... and so on, so as to obtain the nth row of projection data at the (n-1)*θ angle of the 360-degree side of each transverse slice. For another example, taking θ=0.5° and N=720 as an example, a Radon transform is performed starting from an angle of 0° on the 360-degree side of each transverse slice to obtain the first row of projection data, and then a Radon transform is performed at an angle of 0.5° on the 360-degree side of each transverse slice to obtain the second row of projection data, and then a Radon transform is performed at an angle of 1° on the 360-degree side of each transverse slice to obtain the third row of projection data... and so on, so as to obtain the nth row of projection data at an angle of (n-1)*θ on the 360-degree side of each transverse slice.

[0115] It can be understood that a 360-degree Radon transform is performed on all transverse slices of the three-dimensional model. The Radon transform is the projection of the intensity of each transverse slice binary image (two-dimensional grayscale image) along a radial line at a specific angle. Specifically, the Radon transform of each transverse slice binary image is the sum of the Radon transforms of each pixel therein. For example, combined with Figure 3 As shown, the Radon transform operation process for each transverse slice binary image is as follows:

[0116] First, each pixel of the horizontal slice binary image is divided into four sub-pixels, and each sub-pixel is projected separately;

[0117] The contribution of each subpixel is split proportionally to the two nearest bins, based on the distance between the projection position and the bin center;

[0118] Calculate the sum of pixel values ​​based on the projection of the sub-pixel to the center of the bin;

[0119] The situation where the sub-pixel is projected to the center point of the bin is specifically:

[0120] (1) When a sub-pixel is projected to the center of a bin, the bin on the axis will obtain the full value of the sub-pixel, which is one-quarter of the pixel value;

[0121] (2) When a sub-pixel is projected onto the boundary between two bins, the sub-pixel value is evenly split between the two bins.

[0122] Understandably, introducing Figure 3 The Radon transform for a two-dimensional grayscale image shown is already familiar to those skilled in the art and will not be described in detail herein.

[0123] Specifically, in the above step S3, the projection data of the same angle of all the transverse slices of the three-dimensional model are stacked in sequence along the Z-axis direction to form projection slices of the same angle, thereby obtaining projection slices of the side of the three-dimensional model at different angles of 360 degrees. Figure 5a As shown, Figure 5a The generation (transformation) process of projection slices is shown. Specifically, the three-dimensional object is sliced ​​horizontally to obtain multiple layers (sheets) of horizontal slices, and then each layer (sheet) of horizontal slices is subjected to Laden transformation to obtain projection data, and then all horizontal slices are transformed at the same angle (for example, Figure 5a The projection data (preferably after filtering) of the same angle (for example, Figure 5a Projected slices at 0° are shown.

[0124] In this embodiment, different parameters β can be subjected to the printing effect evaluation operation of the embodiment of the present invention to obtain an optimal window function, and filtered back projection transformation can be performed based on the optimal window function.

[0125] Next, how to obtain the optimal window function through the printing effect evaluation operation of the embodiment of the invention will be described in detail. Specifically, in step S4, the optimal window function is determined by the following steps:

[0126] S41, taking M1 transverse slices obtained by equidistant transverse slicing of the three-dimensional model along the Z axis as transverse slice samples; wherein each transverse slice sample is a transverse slice binary image, and M1≥1;

[0127] S42, performing a 360-degree Laden transformation on each of the transverse slice samples, so that each of the transverse slice samples obtains W lines of projection data; wherein each line of projection data is projection data of an angle of each of the transverse slice samples obtained by Laden transformation, and among the W lines of projection data obtained by performing Laden transformation on the same transverse slice sample, the wth line of projection data is projection data of an angle of (w-1)*θ of the transverse slice sample obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°;

[0128] S43, performing filtering back-projection transformation on the W lines of projection data obtained from each transverse slice sample using a to-be-determined window function, thereby obtaining M1 back-projection reconstruction images with the same resolution as each transverse slice sample; wherein the back-projection reconstruction image is the cumulative distribution of light intensity of the transverse slice at the corresponding height of the printed object in the volume bioprinting projection;

[0129] S44, normalizing the pixel value of each back-projection reconstruction image to an integer between 0 and 255 to obtain M1 normalized reconstruction images;

[0130] S45, setting a variable parameter critical curing brightness value, respectively setting the critical curing brightness value to different integers of [1,254], and setting the uncured pixels in each normalized reconstructed image to 0 and the cured pixels to 1 according to the critical curing brightness value, thereby correspondingly obtaining M1*254 simulated printed transverse slices; wherein, the pixels in the normalized reconstructed image that are less than the critical curing brightness value are uncured pixels, and the pixels that are greater than or equal to the critical curing brightness value are cured pixels;

[0131] S46, calculating the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height by the following formula (1), to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height:

[0132]

[0133] Where P is the similarity evaluation value, X and Y are the horizontal resolution and vertical resolution of the horizontal slice sample and the simulated printed horizontal slice; f o (x, y) is the value of the pixel at the position (x, y) in the horizontal slice sample with the upper left corner as the coordinate (1, 1) and the positive direction as downward and left. f(x, y) is the value of the pixel at the position (x, y) in the horizontal slice of the simulated print with the upper left corner as the coordinate (1, 1) and the positive direction as downward and left.

[0134] S47, respectively calculating the average values ​​of M similarity evaluation values ​​corresponding to M1 simulated printed transverse slices at the same critical curing brightness value, and taking the average value closest to 0 as the printing evaluation index of the to-be-determined window function;

[0135] S48, using different window functions to perform filtered back projection transformation on the W lines of projection data obtained for each of the transverse slice samples in step S43, and comparing the printing evaluation indicators of the different window functions obtained after steps S44 to S47, and taking the window function corresponding to the printing evaluation indicator closest to 0 among the printing evaluation indicators of different window functions as the optimal window function.

[0136] First, in step S41, step S41 is a sampling step, that is, a number of transverse slices (for example, M1 slices) are taken out from the transverse slices obtained by equidistant transverse slicing of the three-dimensional model along the Z axis as transverse slice samples, and the printing effect evaluation method of the embodiment of the present invention is performed on the M1 transverse slice samples to obtain an evaluation result.

[0137] It can be understood that in this embodiment, the M1 transverse slice samples used as sampling are Figure 2 The bone screw model shown is divided into 500 slices of equal thickness, and transverse slices corresponding to the heights of the 50th, 150th, 250th, 350th and 450th layers are sampled.

[0138] It can be understood that in step S42, each angle of the 360-degree side of each transverse slice sample is sequentially subjected to Radon transformation, so as to obtain the projection data of the Radon transformation of the corresponding angle of each transverse slice sample. All transverse slice samples (for example, M1 sheets) are subjected to 360-degree Radon transformation, and Radon transformation is the projection of the intensity of each transverse slice sample binary image (two-dimensional grayscale image) along the radial line of a specific angle.

[0139] In this embodiment, after the projection data of each angle of each transverse slice sample is obtained in step S42, a filter back-projection operation is performed on the projection data of each angle of each transverse slice in step S43, so as to obtain a back-projection reconstruction image with the same resolution as that of each transverse slice sample. It can be understood that the back-projection reconstruction image is the cumulative distribution of light intensity of the transverse slice at the corresponding height (e.g., the 150th layer) of the printed object in the volume bioprinting projection.

[0140] Returning to the embodiment of the present invention, continue to combine reference Figure 4 , Figure 7 to Figure 9 Preferably, in step S43, the W lines of projection data obtained from each transverse slice sample are subjected to Fourier transform, filtering, and inverse Fourier transform operations in the filtered back projection transform using a window function to be determined, specifically including:

[0141] S431, performing a fast Fourier transform on the W lines of projection data obtained from each of the transverse slice samples; wherein the width N used in the fast Fourier transform is greater than or equal to the number of pixels of each line of projection data obtained from each of the transverse slice samples;

[0142] S432, using the to-be-determined window function to perform a windowing operation on the ramp filter, and using the ramp filter after the windowing operation to filter each line of projection data of each transverse slice sample after fast Fourier transform;

[0143] S433, performing inverse fast Fourier transform on each line of projection data of each transverse slice sample after filtering, so as to obtain M1 back-projection reconstruction images with the same resolution as that of each transverse slice sample.

[0144] Wherein, the step S432 specifically includes:

[0145] S4321, discretize sampling the ramp filter in the frequency domain; the number of sampling points is an even number greater than or equal to N and closest to N;

[0146] S4322, performing point-to-point multiplication of the discretized ramp filter and the undetermined window function with the same number of discretized sampling points to perform a windowing operation;

[0147] S4323, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head;

[0148] S4324, perform point-to-point multiplication of the ramp filter after windowing and translation operations with each line of projection data after each of the transverse slice samples undergoing fast Fourier transform to perform a filtering operation, and discard the excess points between the ramp filter after windowing and translation operations and the points of each line of projection data after each of the transverse slice samples undergoing fast Fourier transform.

[0149] It can be understood that the use of the above-mentioned filter back-projection transform to filter the projection data of each angle of each transverse slice can effectively reduce star artifacts. Specifically, the Fourier transform of the projection data at each angle is a straight line passing through the center of the frequency domain coordinate, and finally forms a point scattering shape. The density of the center segment at the origin of the ωx-ωy plane is higher than the density in the area far away from the origin, and the area near the origin of the Fourier space is a low-frequency area. Excessive weighting of low-frequency components causes star artifacts in the image. In order to reduce star artifacts, the Fourier space is weighted and corrected to make its density uniform. Therefore, a low-frequency filter, |w| (i.e., a ramp filter) is used to suppress low-frequency components, reduce star artifacts, and improve image clarity.

[0150] In addition, in step S4323, it is necessary to perform a translation operation on the windowed ramp filter data, taking into account that the data after the fast Fourier transform and the windowed ramp filter data do not correspond one to one in order, because the data after the fast Fourier transform is the first half (when the width N is an odd number, it is the first (N+1) / 2 data) representing the data from low frequency to high frequency, the subsequent data is mirror-symmetrical to the previous data.

[0151] It can be understood that in step S43, as a preferred solution, according to the basic principle of filter back-projection, only the projection data with a projection data angle range greater than or equal to 180 degrees in the W lines of projection data obtained from each of the transverse slice samples (that is, there is no need to use the projection data at all angles) can be subjected to filter back-projection transformation using a yet-to-be-determined window function, thereby obtaining M1 back-projection reconstruction images with the same resolution as that of each transverse slice sample.

[0152] Further, in step S43, the undetermined window function refers to a plurality of different window functions used to participate in the evaluation of the printing effect evaluation method provided by the embodiment of the present invention to finally determine the optimal window function. The different window functions refer to Kaiser window functions under different parameters β, and the Kaiser window function formula is shown in formula (3):

[0153]

[0154] Where n = 1, 2, 3, ..., N-1, N represents the total length of the window function; I 0 represents the first kind of Bessel function; β is a variable parameter. For example, in this embodiment, the parameter β can be set to 11 different parameters with values ​​of 0 to 10 to obtain corresponding 11 different window functions for evaluation.

[0155] Furthermore, in step S44 to step S45, in order to obtain the simulated printed transverse slice of the "corresponding height" of the printed object, the pixel values ​​of the back-projected reconstruction image are first normalized to integers from 0 to 255 to obtain a normalized reconstruction image. Then, a variable parameter "critical curing brightness value" is set, ranging from any integer from 1 to 254. The critical curing brightness value represents: the part below the critical curing brightness value in the normalized reconstruction image belongs to the uncured part, that is, the part without the printed object, and the part above or equal to the critical curing brightness value belongs to the cured part, that is, the part with the printed object formed. Through the critical curing brightness value, the uncured pixels in the normalized reconstruction image are set to 0, and the cured pixels are set to 1, so as to obtain a binary image of the transverse slice of the "corresponding height" of the printed object (i.e., the simulated printed transverse slice). It can be seen from the volumetric bioprinting projection process that the transverse slice is a binary image, in which the pixel value of 1 represents the three-dimensional model, and the pixel value of 0 represents the empty.

[0156] refer to Figure 4 , Figure 4 The figure shows the transverse slice samples (original binary images) of the 50th and 200th layers of the three-dimensional model (bone screw model), the back projection images (filtered back projection grayscale images) corresponding to the 50th and 200th transverse slice samples obtained by executing step 3 using the Kaiser window function with parameter β=5, and the final printing effect images corresponding to the 50th and 200th transverse slices obtained by executing step S7 when the critical curing brightness value is set to 115.

[0157] Next, combine Figure 7 to Figure 9 , steps S46 to S48 of the embodiment of the present invention are described in detail. Figure 7 This is a line graph showing the average evaluation value of the Kaiser window function with different β values ​​and the change in critical curing light intensity. Figure 7 In the step S4, the optimal window function is obtained by sampling transverse slices of 50, 150, 250, 350, and 450 layers of bone screws, and the Kaiser window function with different β values ​​and the average evaluation value (printing evaluation index) of the sampling transverse slice evaluation values ​​under different critical curing light intensity values ​​are obtained. Figure 8 This is a line graph of the optimal evaluation value of the Kaiser window function with different β values. Figure 7 ,exist Figure 8 For Kaiser window functions with different β values, the smallest average evaluation value of 254 different critical curing light intensity values ​​is taken as the optimal evaluation value of the Kaiser window function with this β value. Figure 8 It can be seen that when β = 0, 1, 2, 3, 4, and 5, the optimal evaluation value is the smallest 0. Here, the Kaiser window function with β = 5 is taken as the optimal window function. Fig. 9This is a line graph of the average evaluation value of the Kaiser window function with β=5 versus the critical curing light intensity value. When β=5, the average evaluation value of different critical curing light intensity values. Here, the part with the smallest average evaluation value (critical curing light intensity values ​​105-151) is intercepted. It can be seen that when the critical curing light intensity values ​​are 115, 116, and 117, the average evaluation value is the smallest 0. The embodiment of the present invention uses 115 as the critical curing light intensity value.

[0158] Specifically, in steps S46 to S47, the simulated printed transverse slices are calculated with the transverse slice samples of the corresponding height by the above formula (1), and the similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height are obtained. When performing this operation, the transverse slice samples at the same height are compared with the pixels at the same position of the simulated printed transverse slices. If they are different, 1 is added, and then the ratio of the final result to the total number of pixels (the total number of pixels of the transverse slice samples or the simulated printed transverse slices) is calculated. It can be seen from formula (1) that the closer the similarity evaluation value P is to 0, the higher the similarity between the transverse slice samples at the same height and the simulated printed transverse slices, and vice versa.

[0159] As another optional implementation scheme, in step S46, the following formula (2) is used instead of formula (1) to calculate the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height:

[0160]

[0161] Among them, when executing the solution of the above formula (2), in the step S3207, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in the step S48, the window function corresponding to the printing evaluation index closest to 1 among the printing evaluation indexes of different parameters is used as the optimal window function.

[0162] For all transverse slice samples (for example, M1 sheets, and each sheet corresponds to a height, for example, the 100th layer), the above formula (1) or formula (2) is used to calculate 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height, and 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height are obtained. Then, the average value of the M1 similarity evaluation values ​​corresponding to all simulated printed transverse slices (for example, M1 simulated printed transverse slices) corresponding to the same critical curing brightness value (the specific value range is 1 to 254) is calculated respectively, and the optimal average value (when formula (1) is used, the closer the average value is to 0, the better; when formula (2) is used, the closer the average value is to 1, the better) is selected as the printing evaluation index of the undetermined window function.

[0163] Then, in step S48, all different window functions that need to be judged as the pending window functions in step S43 are subjected to filtered back projection transformation and after steps S44 to S47, the printing evaluation indicators of the horizontal slices of different window functions are obtained, the printing evaluation indicators of different window functions are compared, and the window function corresponding to the best printing evaluation indicator (when formula (1) is used, the closer the printing evaluation indicator is to 0, the better; when formula (2) is used, the closer the printing evaluation indicator is to 1, the better) is selected as the optimal window function.

[0164] Furthermore, in step S5 of the embodiment of the present invention, the projection slice is filtered based on a filtered back projection algorithm. Specifically, step S5 further includes:

[0165] S51, performing a fast Fourier transform on each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model; wherein the width N used in the fast Fourier transform is greater than or equal to the number of pixels of each line of projection data of each transverse slice;

[0166] S52, using the optimal window function to perform a windowing operation on the ramp filter, and using the ramp filter after the windowing operation to filter each line of projection data of each projection slice after fast Fourier transform;

[0167] S53, performing inverse fast Fourier transform on each line of projection data of each projection slice after filtering, so as to obtain a filtered projection slice.

[0168] Wherein, the step S52 specifically includes:

[0169] S521, discretizing the ramp filter in the frequency domain for sampling; the number of sampling points is an even number greater than or equal to N and closest to N;

[0170] S522, performing point-to-point multiplication of the discretized ramp filter and the optimal window function with the same number of discretized sampling points to perform a windowing operation;

[0171] S523, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head;

[0172] S524, perform point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after the fast Fourier transform of each projection slice to perform a filtering operation, and discard the extra points of the ramp filter after the windowing and translation operation compared with the points of each line of projection data after the fast Fourier transform of each projection slice.

[0173] like Figure 5b As shown, the projection slices obtained by the volume bio-printing light intensity distribution control method based on the projection algorithm provided by the present invention are used to perform volume printing of the three-dimensional model. During the volume printing process of the three-dimensional model, the projection slices are projected sequentially (for example, timed) in angular order into a printing bottle containing the specific bio-ink, and the printing bottle is rotated at a constant speed so that the angle of the projected projection slices corresponds to the angle of the printing bottle.

[0174] Specifically, the projection slice is projected from a certain angle into a printing bottle containing the specific biological ink in a clockwise or counterclockwise order, and the printing bottle is rotated at a constant speed so that the angle of the projected projection slice corresponds to the angle of the printing bottle. The projection light passes through the printing bottle at a corresponding angle and performs a similar back-projection effect on all transverse sections in each printing bottle. That is, a string of projection data at each angle (here, the grayscale of each row of pixels) is "wiped back" along the original angle and averaged over all two-dimensional pixel points, which will cause the cross-section light intensity in the printing bottle to accumulate into the shape of the corresponding three-dimensional object cross-section, thereby solidifying this cross-section into a three-dimensional object cross-section slice. Light intensity accumulation refers to the superposition of light intensity energy at the same position during the entire printing process. Because the projection slices at each angle are formed by stacking together the projection data of all the cross-sections of the three-dimensional object after filtering at this angle (the projection data of the cross-section of the three-dimensional object after filtering at this angle is a string of data, and this string of data is arranged row by row from low to high in the order of the corresponding cross-sections from bottom to top), so when it is projected into the printing bottle, a similar back-projection operation will be performed on all cross-sections at the same time. The bio-ink of all cross-sections in the printing bottle is simultaneously photocured into the corresponding cross-sections of the three-dimensional object, which is the three-dimensional object from a three-dimensional perspective, thus completing the volume printing of the 3D model.

[0175] Next, combine Figure 6 and Fig.10 Steps S31 to S33 are described in detail. Figure 6 yes Figure 2 A schematic diagram of a 0 degree unfiltered projection slice of the bone screw 3D model is shown. Fig.10 yes Figure 2 The bone screw 3D model shown is a schematic diagram of a 0-degree filtered projection slice after being filtered using a ramp filter windowed with a Kaser window function of β=5.

[0176] Specifically, the volume bioprinting light intensity distribution control method based on the projection algorithm provided by the embodiment of the present invention is to add the Fourier transform, filtering and inverse Fourier transform in the filter back-projection algorithm to the volume bioprinting, and use the above-mentioned filter back-projection transform to filter each line of projection data of the projection slice at different angles of 360 degrees on the side of the three-dimensional model, which can effectively reduce the star-shaped artifacts. Specifically, the Fourier transform of each line of projection data of each angle projection slice is a straight line passing through the center of the frequency domain coordinate, and finally forms a point scattering shape. The density of the center segment at the origin of the ωx-ωy plane is higher than the density in the area far away from the origin, and the area near the origin of the Fourier space is a low-frequency area. Excessive weighting of low-frequency components causes star-shaped artifacts in the image. In order to reduce star-shaped artifacts, the Fourier space is weighted and corrected to make its density uniform. Therefore, a low-frequency filter, |w| (i.e., a ramp filter) is used to suppress low-frequency components, reduce star-shaped artifacts, and improve image clarity.

[0177] In addition, in step S523, it is necessary to perform a translation operation on the windowed ramp filter data, taking into account that the data after the fast Fourier transform and the windowed ramp filter data do not correspond one to one in order, because the data after the fast Fourier transform is the first half (when the width N is an odd number, it is the first (N+1) / 2 data) representing the data from low frequency to high frequency, the subsequent data is mirror-symmetrical to the previous data.

[0178] in, Fig.10 The 0 degree filtered slice after being filtered by the ramp filter using the Kaiser window function with parameter β=5 is shown. It can be understood that the Kaiser window function with different parameters β is shown in formula (3):

[0179]

[0180] Where n = 1, 2, 3, ..., N-1, N represents the total length of the window function; I 0 represents the Bessel function of the first kind; β is a variable parameter.

[0181] Correspondingly, the formula of the Kaiser window function with parameter β = 5 is as follows:

[0182]

[0183] Further, step S6 of the volume bio-printing light intensity distribution control method based on projection algorithm provided in an embodiment of the present invention is described in detail. It can be understood that step S6 is to obtain the filtered projection slice and adjust the pixel value of the filtered projection slice before performing volume printing, so as to make the specific bio-ink used in the printing process solidify when the pixel value of the projection slice is greater than or equal to the critical solidification brightness value.

[0184] The specific execution process of the pixel value adjustment algorithm in step S62 is as follows:

[0185] Multiplying the pixel values ​​of all pixels in the scaled projection slice by a pixel value adjustment parameter (i.e., a printing light intensity parameter), and rounding the result to an integer;

[0186] The pixel values ​​greater than 255 in the obtained results are correspondingly modified to 255, and the modified results are used as new pixel values ​​of the corresponding pixels, so as to obtain a projection slice with adjusted pixel values, and the projection slice with adjusted pixel values ​​is used as a pre-projection slice.

[0187] The following are some examples:

[0188] When printing a specific volume, the bio-ink used in the printing process has a critical light intensity value (collectively referred to as the critical curing light intensity value in this embodiment) when it is fixed. This embodiment performs a printing light intensity adjustment operation on the filtered projection slice, the purpose of which is to use the critical curing brightness value (i.e., pixel value, specifically grayscale value in this embodiment) corresponding to the above-mentioned optimal window function. For example, the critical curing brightness value in this embodiment is 115, such as Figure 4 The light intensity projected by the pixel of the filtered projection slice (as shown) is adjusted to the critical curing light intensity value of the biological ink. There are two adjustment methods, one is to adjust the pixel value (grayscale value) of the filtered projection slice as a whole, and the other is to adjust the light intensity of the projector. In this embodiment, the printing light intensity adjustment operation is performed in combination with the two methods. For example, if in 3D volume bioprinting, the light intensity currently projected by the projector at a grayscale value of 140 is the critical curing light intensity value of the biological ink used, first adjust the light intensity of the projector so that the light intensity projected by the projector at an adjusted pixel value (for example, 118) close to the critical curing brightness value (115) is the critical curing light intensity value; then, adjust the grayscale values ​​of all pixel points of the filtered projection slice through the above-mentioned pixel value adjustment algorithm, and adjust (increase or decrease) the grayscale value as a whole so that the grayscale value of the pixel point corresponding to the critical curing brightness value (115) is adjusted (increased here) to the adjusted pixel value 118.

[0189] It is understandable that before obtaining the filtered projection slice and performing the printing light intensity adjustment operation in step S6, the size of the filtered projection slice needs to be adjusted first, so that the size-adjusted projection slice matches the size of the projection screen used in the printing process. Specifically, Fig.11 As shown, Fig.11 yes Figure 2 The bone screw 3D model shown is a schematic diagram of a scaled projection slice after filtering at 0 degrees and removing invalid parts after using a ramp filter with a Kaser window function of β=5. The printing size adjustment process specifically includes the following steps:

[0190] S6011, performing a frame removal operation on each of the filtered projection slices to obtain a valid projection slice; the frame removal operation includes removing all black columns on the left and right sides, removing all black rows on the top and bottom, and retaining only a valid area in the middle;

[0191] S6012, according to the magnification or reduction factor, multiply the magnification or reduction factor by the side length of the effective projection slice to obtain a scaled projection slice size;

[0192] S6013. Use an image scaling algorithm (for example, a bi-triple interpolation algorithm) to scale the valid projection slice according to the magnification or reduction factor. When zooming in, if the size of the scaled projection slice exceeds the size of the projection screen, calculate the portion of the middle valid area that is the same size as the projection screen, and use the calculated projection slice as the scaled projection slice.

[0193] After the filtered projection slices are preprocessed by adjusting the printing size and the printing light intensity to obtain pre-projection slices, a projection control process is then performed to perform volume printing on the three-dimensional model based on the pre-projection slices. Figure 5b The projection control process of step S7 specifically includes:

[0194] S71, start the stepper motor, control the stepper motor to rotate at a constant speed of K degrees per second, so as to drive the printing bottle arranged on the stepper motor to rotate synchronously; wherein the printing bottle contains the specific biological ink; wherein 0<K, preferably K=6;

[0195] S72, start the projection device to project the pre-projected slice onto the side of the printing bottle, so that the position of the pre-projected slice on the projection screen satisfies:

[0196] X=[(x 0 -x)÷2]

[0197] Y=[(y 0 -y)÷2]

[0198] Wherein, (X, Y) is the coordinate of the upper left corner of the pre-projected slice on the projection screen, and the coordinate is calculated by taking the upper left corner of the projection screen as the origin (0, 0), the horizontal rightward direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; 0 ,y 0 are the width and height of the projection screen; x, y are the width and height of the pre-projection slice; it can be understood that the position coordinates (X, Y) of the upper left corner of the pre-projection slice on the projection screen can be calculated using the above formula, so that the pre-projection slice is always located in the center of the projection screen; wherein, when the projection device is started to start projection, the timer is started at the same time to start timing.

[0199] S73, controlling the pre-projected slices to be projected onto the side of the printing bottle in angular order, and the projection speed is the same as that of the stepping motor, so that the angle between the pre-projected slices and the printing bottle always remains the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

[0200] S74, adjusting the brightness of the projected three-color light (visible light) by using the red, green and blue light intensity parameters;

[0201] S75. When the timer reaches the preset printing time parameter, printing is stopped; wherein, when printing is stopped, the projection is first turned off, and then the stepper motor is stopped.

[0202] The above disclosures are only some preferred embodiments of the present invention, which certainly cannot be used to limit the scope of rights of the present invention. Ordinary technicians in this field can understand that all or part of the processes of the above embodiments and equivalent changes made according to the claims of the present invention are still within the scope of the invention.

Claims

1. A method for controlling light intensity distribution in volumetric bioprinting based on projection algorithm, characterized in that: The method comprises the steps of: S1, slicing the three-dimensional model horizontally along the Z axis at equal intervals to obtain M horizontal slices; wherein each horizontal slice obtained is a horizontal slice binary image, M≥1; S2, performing a 360-degree Laden transformation on each of the transverse slices, so that each transverse slice obtains W lines of projection data; wherein each line of projection data is projection data of an angle of each transverse slice obtained by Laden transformation, and among the W lines of projection data obtained by performing Laden transformation on the same transverse slice, the wth line of projection data is projection data of an angle of (w-1)*θ of the transverse slice obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; S3, stacking the projection data of the same angle of the M transverse slices in sequence along the direction of the Z axis to form a projection slice of the same angle, thereby obtaining projection slices of the side of the three-dimensional model at different angles of 360 degrees; S4, taking any M1 of the M transverse slices as transverse slice samples and executing a volumetric bioprinting evaluation algorithm to obtain an optimal window function; wherein the optimal window function has a corresponding critical curing brightness value; S5, using the optimal window function and based on the filtered back projection algorithm, filtering each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model, so as to obtain filtered projection slices; S6, performing a printing light intensity adjustment operation on the filtered projection slice, by simultaneously adjusting the printer light intensity and adjusting the pixel value of the filtered projection slice so that the specific bio-ink is cured when the pixel value of the projection slice is greater than or equal to the critical curing brightness value, thereby obtaining a pre-projection slice; Wherein, the step S6 comprises: S61, adjusting the light intensity of the projector of the printer so that the light intensity projected by the projector when the pixel value is close to the critical curing brightness value is equal to the critical curing light intensity value of the specific biological ink; S62, adjusting the pixel values ​​of all pixel points of the filtered projection slice by the following pixel value adjustment algorithm, so that the pixel value of the pixel point corresponding to the critical solidification brightness value is adjusted to the adjusted pixel value, and using the projection slice after the pixel value adjustment as the pre-projection slice: ； Among them, h is the adjusted pixel value, b is the pixel value adjustment parameter, is the original pixel value, <> represents rounding; S7, applying the pre-projected slices to perform volumetric bio-printing, projecting the pre-projected slices into a printing bottle containing the specific bio-ink in order of angle, and rotating the printing bottle at a constant speed so that the angle of the projected slices corresponds to the angle of the printing bottle to perform volumetric printing of the three-dimensional model; wherein the bio-inks of all transverse sections in the printing bottle are simultaneously photocured into corresponding transverse slices of the three-dimensional object; Wherein, the step S4 specifically includes: S41, taking M1 transverse slices obtained by equidistant transverse slicing of the three-dimensional model along the Z axis as transverse slice samples; wherein each transverse slice sample is a transverse slice binary image, and M1≥1; S42, performing a 360-degree Laden transformation on each of the transverse slice samples, so that each of the transverse slice samples obtains W lines of projection data; wherein each line of projection data is projection data of an angle of each of the transverse slice samples obtained by Laden transformation, and among the W lines of projection data obtained by performing Laden transformation on the same transverse slice sample, the wth line of projection data is projection data of an angle of (w-1)*θ of the transverse slice sample obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; S43, performing filtering back-projection transformation on the W lines of projection data obtained from each transverse slice sample using a to-be-determined window function, thereby obtaining M1 back-projection reconstruction images with the same resolution as each transverse slice sample; wherein the back-projection reconstruction image is the cumulative distribution of light intensity of the transverse slice at the corresponding height of the printed object in the volume bioprinting projection; S44, normalizing the pixel value of each back-projection reconstructed image to an integer between 0 and 255 to obtain M1 normalized reconstructed images; S45, setting a variable parameter critical curing brightness value, respectively setting the critical curing brightness value to different integers of [1, 254], and setting the uncured pixels in each normalized reconstructed image to 0 and the cured pixels to 1 according to the critical curing brightness value, thereby correspondingly obtaining M1*254 simulated printed transverse slices; wherein, in the normalized reconstructed image, pixels less than the critical curing brightness value are uncured pixels, and pixels greater than or equal to the critical curing brightness value are cured pixels; S46, calculating the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height by the following formula (1), to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height: Formula (1): Among them, P is the similarity evaluation value, X and Y are the horizontal resolution and vertical resolution of the horizontal slice sample and the simulated printed horizontal slice; It is the value of the pixel at the position (x,y) in the horizontal slice sample with the upper left corner as the coordinate (1,1) and the positive direction as downward and left; It is the value of the pixel at the position (x,y) in the simulated printed horizontal slice with the upper left corner as the coordinate (1,1) and the positive direction downward and to the left as the reference; S47, respectively calculating the average values ​​of M1 similarity evaluation values ​​corresponding to M1 simulated printed transverse slices at the same critical curing brightness value, and taking the average value closest to 0 as the printing evaluation index of the to-be-determined window function; S48, using different window functions to perform filtered back projection transformation on the W lines of projection data obtained for each of the transverse slice samples in step S43, and comparing the printing evaluation indicators of the different window functions obtained after steps S44 to S47, and taking the window function corresponding to the printing evaluation indicator closest to 0 among the printing evaluation indicators of different window functions as the optimal window function.

2. According to the projection algorithm-based volume bioprinting light intensity distribution control method of claim 1, it is characterized in that: The step S7 comprises: S71, start the stepper motor, control the stepper motor to rotate at a constant speed of K degrees per second, so as to drive the printing bottle arranged on the stepper motor to rotate synchronously; wherein the printing bottle contains the specific biological ink, 0 <K; S72, start the projection device to project the pre-projected slice onto the side of the printing bottle, so that the position of the pre-projected slice on the projection screen satisfies: ； ； Wherein, (X, Y) is the coordinate of the upper left corner of the pre-projected slice on the projection screen, and the coordinate is calculated by taking the upper left corner of the projection screen as the origin (0, 0), the horizontal rightward direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; , is the width and height of the projection screen; , is the width and height of the pre-projected slice; S73, controlling the pre-projected slices to be projected onto the side of the printing bottle in angular order, and the projection speed is the same as that of the stepping motor, so that the angle between the pre-projected slices and the printing bottle always remains the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

3. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 1, characterized in that: A print size adjustment step is also included between step S5 and step S6, and the print size adjustment step includes: Performing a border removal operation on each of the filtered projection slices to obtain a valid projection slice; the border removal operation includes removing all black columns on the left and right sides, removing all black rows on the top and bottom sides, and retaining only the valid area in the middle; According to the magnification or reduction factor, multiplying the magnification or reduction factor by the side length of the effective projection slice to obtain the scaled projection slice size; Using an image scaling algorithm to scale the effective projection slice according to the magnification or reduction multiple, when scaling up, if the size of the scaled projection slice exceeds the size of the projection screen, then calculating the portion of the middle effective area that has the same size as the projection screen, and using the calculated projection slice as the scaled projection slice; The filtered projection slice in step S6 is the scaled projection slice.

4. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 3, characterized in that: The image scaling algorithm is a bi-triple interpolation algorithm.

5. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 2, characterized in that: In the step S7, a printing time control step is also included, wherein the printing time control step starts timing when the projection device is started to start projection, and stops printing when the printing time parameter is reached; wherein, when stopping printing, the projection is first turned off, and then the rotation of the stepper motor is stopped.

6. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 1, characterized in that: In step S46, the following formula (2) is used instead of formula (1) to calculate the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height, so as to obtain 254 similarity evaluation values ​​corresponding to the 254 simulated printed transverse slices at each height: Formula (2): Wherein, in said S47, the average value closest to 1 is taken as the printing evaluation index of the undetermined window function; in said S48, the window function corresponding to the printing evaluation index closest to 1 among the printing evaluation indexes with different parameters is taken as the optimal window function.

7. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 1, characterized in that: The step S5 specifically includes: S51, performing a fast Fourier transform on each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model; wherein the width N used in the fast Fourier transform is greater than or equal to the number of pixels of each line of projection data of each transverse slice; S52, using the optimal window function to perform a windowing operation on the ramp filter, and using the ramp filter after the windowing operation Filtering each line of projection data of each projection slice after fast Fourier transformation; S53, performing inverse fast Fourier transform on each line of projection data after filtering of each projection slice, The filtered projection slice is obtained.

8. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 7, characterized in that: The step S52 specifically includes: S521, discretizing the ramp filter in the frequency domain for sampling; the number of sampling points is an even number greater than or equal to N and closest to N; S522, performing point-to-point multiplication of the discretized ramp filter and the optimal window function with the same number of discretized sampling points to perform a windowing operation; S523, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head; S524, perform point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after the fast Fourier transform of each projection slice to perform a filtering operation, and discard the extra points of the ramp filter after the windowing and translation operation compared with the points of each line of projection data after the fast Fourier transform of each projection slice.

9. The method for controlling light intensity distribution of volumetric bioprinting based on projection algorithm according to claim 1, characterized in that: In step S2, Laden transformation is performed on each angle of the 360-degree side of each transverse slice in sequence, so as to obtain Laden transformed projection data of the corresponding angle of each transverse slice.

Images