A volumetric bioprinting control method

Abstract

The present invention provides a volume bioprinting control method, comprising: slicing a three-dimensional model in equidistant transverse directions along the Z axis to obtain M transverse slices; performing a 360-degree Laden transform on each transverse slice to obtain W lines of projection data for each transverse slice; performing a filtering operation on each line of projection data, and stacking the 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, thereby obtaining projection slices at different angles of 360 degrees on the side of the three-dimensional model; projecting the obtained projection slices in sequence according to the angle sequence into a printing bottle filled with photocurable bio-ink, and rotating the printing bottle at a uniform speed so that the angle of the projected projection slice corresponds to the angle of the printing bottle to perform volume 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.

Description

Technical Field

[0001] The present invention belongs to the technical field of 3D printing, and in particular relates to a volume bioprinting control method. 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 control method, device, equipment and storage medium, 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 volume bioprinting control method, 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, performing a filtering operation on each line of projection data, and stacking the 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, so as to obtain projection slices at different angles of 360 degrees on the side of the three-dimensional model;

[0010] S4, projecting the obtained projection slices into a printing bottle containing photocurable biological ink in order of angles, and rotating the printing bottle at a constant speed so that the angle of the projected projection slices corresponds to the angle of the printing bottle to perform volume printing of the three-dimensional model; wherein the biological ink of all transverse sections in the printing bottle is simultaneously photocured into corresponding transverse slices of the three-dimensional object;

[0011] The step S4 comprises:

[0012] S41, starting a stepper motor, controlling the stepper motor to rotate at a constant speed of K degrees per second, thereby driving a printing bottle disposed on the stepper motor to rotate synchronously; wherein the printing bottle contains a photocurable biological ink;

[0013] S42, start the projection device to project the projection slice onto the side of the printing bottle, so that the position of the projection slice on the projection screen satisfies:

[0014] X=[(x0-x)÷2]

[0015] Y=[(y0-y)÷2]

[0016] Wherein (X, Y) is the coordinate of the upper left corner of the projection 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 right direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; x0, y0 are the width and height of the projection screen; x, y are the width and height of the projection slice;

[0017] S43, controlling the projection slices to be projected onto the side of the printing bottle in angular order, and the switching speed of the projection slices is the same as that of the stepping motor, so that the angles of the projection slices and the printing bottle are always kept the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

[0018] As an improvement of the above solution, a printing light intensity adjustment step is further included between step S3 and step S4, and the printing light intensity adjustment step includes:

[0019] Multiplying the pixel values ​​of all pixels of the projection slice by the printing light intensity parameter, and rounding the result to an integer;

[0020] Modify the pixel values ​​greater than 255 in the obtained results to 255, and use the modified results as new pixel values ​​of the corresponding pixels, so as to obtain a projection slice with adjusted pixel values, and use the projection slice with adjusted pixel values ​​as a pre-projection slice;

[0021] In step S4, the projected slice is the pre-projected slice; step S4 further comprises:

[0022] S44. Adjust the brightness of the projected visible light through the red, green and blue light intensity parameters.

[0023] As an improvement of the above solution, a printing size adjustment step is further included before the printing light intensity adjustment step, and the printing size adjustment step includes:

[0024] Performing a border removal operation on each of the 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;

[0025] 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;

[0026] 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;

[0027] Wherein, the projection slice in the printing light intensity adjustment step is the scaled projection slice.

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

[0029] As an improvement of the above scheme, step S4 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.

[0030] As an improvement of the above solution, step S3 specifically includes:

[0031] S31, 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;

[0032] S32, using a 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;

[0033] S33, 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.

[0034] As an improvement of the above solution, step S3 specifically includes:

[0035] Wherein, the step S32 specifically includes:

[0036] S321, 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;

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

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

[0039] S324, 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 excess 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.

[0040] As an improvement of the above solution, the window function is an optimal window function, and the optimal window function is determined by the following steps:

[0041] S3201, 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;

[0042] S3202, perform 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°;

[0043] S3203, using a to-be-determined window function to perform filtered back-projection transformation on the W lines of projection data obtained from each of the transverse slice samples, thereby obtaining M1 back-projection reconstruction images with the same resolution as each of the transverse slice samples; 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 volumetric bioprinting projection;

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

[0045] S3205, 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;

[0046] S3206, 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:

[0047]

[0048] 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.

[0049] S3207, respectively calculating the average values ​​of M similarity evaluation values ​​corresponding to M 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;

[0050] S3208. Use 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 S3203, and compare the printing evaluation indicators of the different window functions obtained after steps S3204 to S3207, and use 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.

[0051] As an improvement of the above solution, in step S3206, 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:

[0052]

[0053] Among them, in S3207, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in S3208, the window function corresponding to the printing evaluation index closest to 1 among the printing evaluation indexes with different parameters is used as the optimal window function.

[0054] 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.

[0055] Compared with the prior art, the volume bioprinting control method provided by the embodiment of the present invention has the following technical effects:

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

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

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

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

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

[0061] ③. 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.

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

[0063] ① 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

[0064] 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.

[0065] Figure 1 A schematic flow chart of a volumetric bioprinting control method provided in an embodiment of the present invention.

[0066] Figure 2 A schematic structural diagram of a transverse slice generated by a volume bioprinting control method provided by an embodiment of the present invention.

[0067] Figure 3 A schematic diagram of a Radon transform performed on a transverse slice sample by a volumetric bioprinting control method provided in an embodiment of the present invention.

[0068] Figure 4 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.

[0069] Figure 5 A schematic diagram of volumetric printing of a three-dimensional model of the obtained projection slices using a volumetric bioprinting control method provided by an embodiment of the present invention.

[0070] Figure 6 This is a schematic diagram of a 90-degree unfiltered projection slice of the DNA model.

[0071] Figure 7 It is a schematic diagram of the filtered projection slice of the DNA model after executing steps S3 to S5 at 90 degrees.

[0072] Figure 8The figure shows a comparison between a binary image of a transverse slice corresponding to the 400th transverse slice of the DNA model reconstructed after filtering operation using the prior art and a binary image of a transverse slice reconstructed after filtered back projection using the projection filter of the present invention (ie, a filtered back-projected grayscale image).

[0073] Fig. 9 It is a curve diagram of the change of the printing evaluation index and the critical curing brightness value of the Kaiser window function with different parameters β.

[0074] Fig.10 The printed evaluation index under the Kaiser window function with different parameters β is shown.

[0075] Fig.11 The figure shows the change of the printing evaluation index (average evaluation value) of the Kaiser window function with β being 0 with the critical curing light intensity value.

[0076] Fig.12 Shown is the print sizing process for a 0-degree projection slice of a DNA model. DETAILED DESCRIPTION

[0077] 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.

[0078] Figure 1 A flow chart of a volume bioprinting control method provided in an embodiment of the present application, the method comprising steps S1 to S4:

[0079] 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;

[0080] 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°;

[0081] S3, performing a filtering operation on each line of projection data, and stacking the 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, so as to obtain projection slices at different angles of 360 degrees on the side of the three-dimensional model;

[0082] S4. Projecting the obtained projection slices sequentially into a printing bottle containing photocurable bio-ink in angular order, and rotating the printing bottle at a constant speed so that the angle of the projected projection slices corresponds to the angle of the printing bottle to perform volume printing of the three-dimensional model; wherein the bio-ink of all transverse sections in the printing bottle is simultaneously photocured into corresponding transverse slices of the three-dimensional object.

[0083] Wherein, the step S4 comprises:

[0084] S41, starting a stepper motor, controlling the stepper motor to rotate at a constant speed of K degrees per second, thereby driving a printing bottle disposed on the stepper motor to rotate synchronously; wherein the printing bottle contains a photocurable biological ink;

[0085] S42, start the projection device to project the projection slice onto the side of the printing bottle, so that the position of the projection slice on the projection screen satisfies:

[0086] X=[(x0-x)÷2]

[0087] Y=[(y0-y)÷2]

[0088] Wherein (X, Y) is the coordinate of the upper left corner of the projection 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 right direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; x0, y0 are the width and height of the projection screen; x, y are the width and height of the projection slice;

[0089] S43, controlling the projection slices to be projected onto the side of the printing bottle in angular order, and the switching speed of the projection slices is the same as that of the stepping motor, so that the angles of the projection slices and the printing bottle are always kept the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

[0090] Below, the various steps of a volume bioprinting control method provided in an embodiment of the present application will be described in detail in conjunction with relevant drawings.

[0091] 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.

[0092] 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.

[0093] Specifically, the following is combined Figure 2 , how to obtain the binary image of each horizontal slice of the three-dimensional model will be described in detail. It can be understood that in S1, the three-dimensional model (DNA 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 = z0, where Z represents the coordinate of the Z axis, z0 is a constant, and z0 is initialized to 0; the step S1 specifically includes steps S11 to S15:

[0094] S11. Compare the Z-axis coordinates of the three vertices of each triangle with z0, and select the triangles that satisfy the conditions that there are both vertices with Z coordinates greater than z0 and vertices with Z coordinates less than z0, or there are both vertices with Z coordinates equal to z0.

[0095] S12. For a triangle having both vertices with Z coordinates greater than z0 and vertices with Z coordinates less than z0, when there is a vertex with a Z coordinate equal to z0, connect the remaining two vertices and calculate the intersection of the connected line segment and the plane parallel to the X and Y axes.

[0096] S13. For a triangle without a vertex whose Z coordinate is equal to z0, calculate the intersection of two vertices whose Z coordinates are greater than z0 and the line segments connecting the vertices whose Z coordinates are less than z0 and the plane parallel to the X and Y axes, or calculate the intersection of two vertices whose Z coordinates are less than z0 and the line segments connecting the vertices whose Z coordinates are greater than z0 and the plane parallel to the X and Y axes;

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

[0098]

[0099]

[0100] Among them, the coordinates of point P1 with Z coordinate greater than z0 are (x1, y1, z1), the coordinates of another point P2 with Z coordinate less than z0 are (x2, y2, z2), the intersection point is P, and 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 variables of the X axis, Y axis, and Z axis, and the origin is O. is the vector from O to P, is the vector from O to P1, is the vector from P1 to P2;

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

[0102]

[0103] in, represents the vector from P1 to P;

[0104] Thus, the coordinates of each intersection point P are obtained as Therefore, each triangle obtains one intersection point, one endpoint with a Z coordinate equal to z0, or two intersection points or two endpoints with a Z coordinate equal to z0, and the two points obtained for one triangle are taken as a group.

[0105] S14, connect two points of all point groups to obtain a contour image of a cross section of the three-dimensional model with a Z-axis coordinate of z0 and parallel to the X-axis and Y-axis planes, pixelate the contour image, and make the pixel points inside the contour image white and the pixel points outside the contour image black, thereby obtaining a binary image of the horizontal slice of the current layer. Figure 2 The transverse slice binary images (original binary images) of the 150th transverse slice and the 400th transverse slice of the three-dimensional model are respectively shown.

[0106] 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.

[0107] S15, determine whether the current z0 is less than H, if so, replace the current z0 with the value of z0+h 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.01㎜.

[0108] 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.

[0109] 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.

[0110] 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.

[0111] 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.

[0112] 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:

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

[0114] 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;

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

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

[0117] (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;

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

[0119] 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.

[0120] In this embodiment, after obtaining the projection data of each angle of each transverse slice of the three-dimensional model through step S2, the projection data of each angle of each transverse slice is first filtered, and then the filtered projection data of the same angle of all transverse slices are stacked in sequence along the direction of the Z axis to form projection slices of the same angle, thereby obtaining projection slices at different angles of 360 degrees on the side of the three-dimensional model. It can be understood that by filtering the projection data of each angle of each transverse slice and then stacking them to form projection slices, the image clarity of 3D volume printing can be further improved.

[0121] Specifically, in the above step S3, the filtered projection data of the same angle of all 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 4 As shown, Figure 4 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 4 The projection data (preferably after filtering) of the same angle (for example, Figure 4 Shown are projected slices at 90°).

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

[0123] S31, 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;

[0124] S32, using a 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;

[0125] S33, 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.

[0126] Wherein, the step S32 specifically includes:

[0127] S321, 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;

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

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

[0130] S324, 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 excess 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.

[0131] like Figure 5 As shown, the projection slices obtained by the volume bioprinting control method provided by the present invention are used to perform volume printing of a three-dimensional model. During the volume printing of the three-dimensional model, the projection slices are projected sequentially (for example, timed) in an angular sequence into a printing bottle containing photocurable 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.

[0132] Specifically, the projection slice is projected from a certain angle in a clockwise or counterclockwise order into a printing bottle containing photocurable bio-ink, and the printing bottle is rotated at a constant speed so that the angle of the projected 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 is 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 cross-section of a three-dimensional object. 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.

[0133] Next, combine Figure 6 to Figure 8 Steps S31 to S33 are described in detail. Figure 6 This is a schematic diagram of a 90-degree unfiltered projection slice of the DNA model. Figure 7 It is a schematic diagram of the filtered projection slice of the DNA model after executing steps S31 to S33 at 90 degrees. Figure 8 The figure shows a comparison between a binary image of a transverse slice corresponding to the 400th transverse slice of the DNA model reconstructed after filtering operation using the prior art and a binary image of a transverse slice reconstructed after filtered back projection using the projection filter of the present invention (ie, a filtered back-projected grayscale image).

[0134] Specifically, the volume bioprinting control method provided by the embodiment of the present invention is to add Fourier transform, filtering and inverse Fourier transform in the filter back-projection algorithm to the volume bio 3D printing, 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, weighted correction is performed on the Fourier space 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.

[0135] In addition, in step S323, 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.

[0136] By comparison Figure 6 and Figure 7 It can be seen that Figure 6 What is displayed is a 90-degree projection slice of the DNA model obtained without executing steps S31 to S33, which is 500 pixels high and 800 pixels wide. Figure 7 Will Figure 6 The 90-degree unfiltered projection slice of the DNA model shown is a 90-degree filtered projection slice of the DNA model after the filtering of steps S31 to S3 by the Kaiser window function with β=0. Figure 7 It can be seen that the edges are clear and the star-shaped artifacts are effectively reduced.

[0137] Likewise, Figure 8 The 400th layer of the DNA model is shown as a binary image of the transverse slice reconstructed after filtering operation (or only back-projection operation) using the prior art (shown as a back-projection grayscale image in the figure), and it can be clearly seen that the edges are blurred and the star-shaped artifacts are serious. The 400th layer of the transverse slice is reconstructed by back-projecting the projection slice using the projection filtering method provided by the embodiment of the present invention (the filtered back-projection grayscale image in the figure), and it can be seen that the edges are clear and the star-shaped artifacts are effectively reduced.

[0138] in, Figure 8 The Kaiser window function with parameter β=0 is used for filtering back projection transformation. The Kaiser window function with different parameters β is shown in formula (3):

[0139]

[0140] Wherein, n=1, 2, 3, ..., N-1, N represents the total length of the window function; I0 represents the first kind of Bessel function; β is a variable parameter.

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

[0142]

[0143] 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.

[0144] 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, the optimal window function is determined by the following steps:

[0145] S3201, 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;

[0146] S3202, perform 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°;

[0147] S3203, using a to-be-determined window function to perform filtered back-projection transformation on the W lines of projection data obtained from each of the transverse slice samples, thereby obtaining M1 back-projection reconstruction images with the same resolution as each of the transverse slice samples; 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 volumetric bioprinting projection;

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

[0149] S3205, 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;

[0150] S3206, 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:

[0151]

[0152] 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.

[0153] S3207, 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;

[0154] S3208. Use 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 S3203, and compare the printing evaluation indicators of the different window functions obtained after steps S3204 to S3207, and use 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.

[0155] First, in step S3201, step S3201 is a sampling step, that is, taking out a number of horizontal slices (for example, M1 slices) from the horizontal slices obtained by equidistantly slicing the three-dimensional model along the Z axis as horizontal slice samples, and executing the printing effect evaluation method of the embodiment of the present invention on the M horizontal slice samples to obtain an evaluation result.

[0156] It can be understood that in this embodiment, the M1 transverse slice sample used as the sampling is a DNA model (see Figure 2 ) is divided into 500 slices of equal thickness, and the transverse slices corresponding to the heights of the 50th, 100th, 150th, 100th, 150th, 200th, 250th, 300th, 350th, 400th and 450th layers are sampled.

[0157] It can be understood that in step S3202, each angle of the 360-degree side of each transverse slice sample is sequentially subjected to Radon transformation, so as to obtain the Radon transformation projection data of the corresponding angle of each transverse slice sample. A 360-degree Radon transformation is performed on all transverse slice samples (e.g., M1 sheets), and Radon transformation is the projection of the intensity of each transverse slice sample binary image (two-dimensional grayscale image) along a radial line of a specific angle.

[0158] In this embodiment, after the projection data of each angle of each transverse slice sample is obtained in step S3202, a filter back-projection operation is performed on the projection data of each angle of each transverse slice in step S3203, 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 100th layer) of the printed object in the volume bioprinting projection.

[0159] Returning to the embodiment of the present invention, continue to refer to Figure 2 , Figure 6 and Figure 7 Preferably, in step S3203, 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:

[0160] S32031, 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;

[0161] S32032, 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;

[0162] S32033. Perform 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.

[0163] Wherein, the step S32032 specifically includes:

[0164] S320321. Discretely sample 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;

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

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

[0167] S320324. Perform point-to-point multiplication of the ramp filter after windowing and translation operations with each line of projection data after fast Fourier transform of each of the horizontal slice samples to perform a filtering operation, and discard the excess points of the ramp filter after windowing and translation operations compared with the points of each line of projection data after fast Fourier transform of each of the horizontal slice samples.

[0168] 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.

[0169] In addition, in step S320323, 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.

[0170] It can be understood that in step S3203, 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.

[0171] Further, in step S3203, 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):

[0172]

[0173] Wherein, n=1, 2, 3, ..., N-1, N represents the total length of the window function; I0 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.

[0174] Furthermore, in step S3204 to step S3205, 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.

[0175] refer to Figure 2 , Figure 2The figures show the 150th and 400th transverse slice samples (original binary images) of the three-dimensional model (DNA model), the back-projection reconstruction images (filtered back-projection grayscale images) corresponding to the 150th and 400th transverse slice samples obtained by executing step 3 using the Kaiser window function (optimal window function) with parameter β=0, and the simulated printed transverse slices (binary images) corresponding to the 150th and 400th transverse slice samples obtained by executing step 5 when the critical curing brightness value is set to 136.

[0176] Next, combine Figures 9 to 11 , the steps S3206 to S3208 of the embodiment of the present invention are described in detail. Fig. 9 As shown, Fig. 9 It is a curve diagram of the change of the printing evaluation index and the critical curing brightness value under the Kaiser window function with different parameters β, showing the change of the printing evaluation index and the critical curing brightness value under the Kaiser window function with different parameters β. Fig. 9 In the figure, the horizontal axis represents different “critical curing brightness values”, and the vertical axis represents the “average similarity evaluation value” of the transverse slice samples, that is, the printing evaluation index. Fig.10 The print evaluation index under the Kaiser window function with different parameters β is shown. Fig.10 In the figure, the horizontal axis represents the Kaiser window function with different β values, and the vertical axis is the printing evaluation index. Fig.11 The figure shows the change of the printing evaluation index (average evaluation value) of the Kaiser window function with β being 0 with the critical curing light intensity value.

[0177] Specifically, in steps S3206 to S3207, 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.

[0178] As another optional implementation scheme, in step S3206, 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:

[0179]

[0180] 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 S3208, 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.

[0181] 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.

[0182] Then, in step S3208, all different window functions that need to be judged as the pending window functions in step S103 are subjected to filtered back projection transformation and after steps S3204 to S3207, 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.

[0183] Next, combine Fig.12 , further describing in detail step S4 of the volume bioprinting control method provided by the embodiment of the present invention. It can be understood that step S4 is to obtain the filtered projection slices, and to perform printing control based on the filtered projection slices to perform volume printing of the three-dimensional model. In specific implementation, in step S4, three processes of printing size adjustment, printing light intensity adjustment and projection control are specifically performed. The following three processes are described in detail in conjunction with the accompanying drawings.

[0184] Before obtaining the filtered projection slices for volume printing, the size of the filtered projection slices needs to be adjusted first, so that the size-adjusted projection slices match the size of the projection screen used in the printing process. Specifically, Fig.12 As shown, Fig.12The print size adjustment process for the 0-degree projection slice of the DNA model is shown. The print size adjustment process specifically includes the following steps:

[0185] S4011, performing a frame removal operation on each of the 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;

[0186] S4012, according to the magnification or reduction factor, multiplying the magnification or reduction factor by the side length of the effective projection slice to obtain a scaled projection slice size;

[0187] S4013. Use an image scaling algorithm (e.g., 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.

[0188] Before obtaining the filtered projection slices for volume printing, after the print size is adjusted, it is necessary to further perform a print intensity adjustment process, which specifically includes the following steps:

[0189] S4021, multiplying the pixel values ​​of all pixels in the scaled projection slice by the printing light intensity parameter, and rounding the result to an integer;

[0190] S4022. Modify the pixel values ​​greater than 255 in the obtained results to 255, and use the modified results as new pixel values ​​of the corresponding pixels, thereby obtaining a projection slice with adjusted pixel values, and use the projection slice with adjusted pixel values ​​as a pre-projection slice.

[0191] It can be understood that the purpose of performing the printing light intensity adjustment process is to make the specific biological 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, specifically by simultaneously adjusting the printer light intensity and adjusting the pixel value of the filtered projection slice so that the specific biological ink solidifies when the pixel value of the projection slice is greater than or equal to the critical solidification brightness value, including:

[0192] 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;

[0193] Adjust the pixel values of all pixel points of the filtered projection slice through the following pixel value adjustment algorithm, adjust the pixel value of the pixel point corresponding to the critical curing brightness value to the adjusted pixel value, and use the projection slice with adjusted pixel values as the pre-projection slice:

[0194]

[0195] where h is the adjusted pixel value, b is the pixel value adjustment parameter (i.e., the printing light intensity parameter), h o is the original pixel value, and <> represents rounding to the nearest integer.

[0196] An example is as follows:

[0197] In specific volume printing, there is a critical light intensity value (collectively referred to as the critical curing light intensity value in this embodiment) when the bio-ink used in the printing process is fixed. The purpose of performing the printing light intensity adjustment operation on the filtered projection slice in this embodiment is to adjust the light intensity projected by the pixels corresponding to the critical curing brightness value (i.e., pixel values, specifically gray values in this embodiment, for example, the critical curing brightness value in this embodiment is 136, as Figure 2 shown) to the critical curing light intensity value of the bio-ink. Among them, there are two adjustment methods. One is to globally adjust the pixel values (gray values) of the filtered projection slice, and the other is to adjust the light intensity of the projector. In this embodiment, the printing light intensity adjustment operation is performed by combining the two methods. For example, if in volume bio-printing, the light intensity projected by the projector at the current 150 gray value is the critical curing light intensity value of the bio-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, 138) close to the critical curing brightness value (136) is the critical curing light intensity value; then, adjust the gray values of all pixel points of the filtered projection slice through the above pixel value adjustment algorithm, and globally adjust (increase or decrease) the gray values so that the gray value of the pixel point corresponding to the critical curing brightness value (136) is adjusted (increased here) to the adjusted pixel value 138.

[0198] After performing preprocessing of printing size adjustment and printing light intensity adjustment on the filtered projection slice to obtain a pre-projection slice, then perform a projection control process to perform volume printing on the three-dimensional model based on the pre-projection slice. The projection control process specifically includes:

[0199] S41. Start the stepping motor, control the stepping motor to rotate at a uniform speed of K degrees per second, so as to drive the printing bottle arranged on the stepping motor to rotate synchronously; wherein, the printing bottle is filled with photocurable bio-ink; wherein, 0 < K, and more preferably K = 6;

[0200] S42, 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:

[0201] X=[(x0-x)÷2]

[0202] Y=[(y0-y)÷2]

[0203] 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 right direction is the positive direction of the X axis, and the vertical downward direction is the positive direction of the Y axis; x0, y0 are the width and height of the projection screen; x, y are the width and height of the pre-projected slice; it can be understood that the above formula can be used to calculate the position coordinate (X, Y) of the upper left corner of the pre-projected slice on the projection screen, so that the pre-projected slice is always located at 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.

[0204] S43, controlling the pre-projected slices to be projected onto the side of the printing bottle in angular order, and the switching speed of the projection slices is the same as that of the stepping motor, so that the angles of the pre-projected slices and the printing bottle are always kept the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

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

[0206] S45, 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 rotation of the stepper motor is stopped.

[0207] 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 volume bioprinting control method, 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, performing a filtering operation on each line of projection data, and stacking the 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, so as to obtain projection slices at different angles of 360 degrees on the side of the three-dimensional model; S4, projecting the obtained projection slices into a printing bottle containing photocurable biological ink in order of angles, and rotating the printing bottle at a constant speed so that the angle of the projected projection slices corresponds to the angle of the printing bottle to perform volume printing of the three-dimensional model; wherein the biological ink of all transverse sections in the printing bottle is simultaneously photocured into corresponding transverse slices of the three-dimensional object; Wherein, the step S3 specifically includes: S31, 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; S32, using a 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; S33, 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; The window function is an optimal window function, and the optimal window function is determined by the following steps: S3201, 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; S3202, perform 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 the transverse slice sample at an angle of (w-1)*θ of the Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; S3203, using a to-be-determined window function to perform filtered back-projection transformation on the W lines of projection data obtained from each of the transverse slice samples, thereby obtaining M1 back-projection reconstruction images with the same resolution as each of the transverse slice samples; 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 volumetric bioprinting projection; S3204, normalizing the pixel value of each back-projection reconstructed image to an integer between 0 and 255, to obtain M1 normalized reconstructed images; S3205, 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; S3206, calculating the 254 simulated printed transverse slices at each height and the transverse slice samples at the corresponding height by the following formula (1), and obtaining 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; S3207, 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; S3208, performing filtering back projection transformation on the W lines of projection data obtained for each of the transverse slice samples in step S3203 using different window functions, and comparing the printing evaluation indicators of the different window functions obtained after steps S3204 to S3207, 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; The step S4 comprises: S41, starting the stepper motor, controlling the stepper motor to rotate at a constant speed of K degrees per second, thereby driving the printing bottle arranged on the stepper motor to rotate synchronously; wherein the printing bottle is filled with photocurable biological ink, 0 <K; S42, start the projection device to project the projection slice onto the side of the printing bottle, so that the position of the projection slice on the projection screen satisfies: ； Wherein, (X, Y) is the coordinate of the upper left corner of the projection 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; , are the width and height of the projection slice; S43, controlling the projection slices to be projected onto the side of the printing bottle in angular order, and the switching speed of the projection slices is the same as that of the stepping motor, so that the angles of the projection slices and the printing bottle are always kept the same; wherein the central axis of the projection screen coincides with the central axis of the printing bottle.

2. A volume bioprinting control method according to claim 1, characterized in that: A printing light intensity adjustment step is also included between step S3 and step S4, and the printing light intensity adjustment step includes: Multiplying the pixel values ​​of all pixels of the projection slice by the printing light intensity parameter, and rounding the result to an integer; Modify the pixel values ​​greater than 255 in the obtained results to 255, and use the modified results as new pixel values ​​of the corresponding pixels, thereby obtaining a projection slice with adjusted pixel values, and use the projection slice with adjusted pixel values ​​as a pre-projection slice; In step S4, the projected slice is the pre-projected slice; step S4 further comprises: S44. Adjust the brightness of the projected visible light through the red, green and blue light intensity parameters.

3. A volume bioprinting control method according to claim 2, characterized in that: Before the printing light intensity adjustment step, a printing size adjustment step is also included, and the printing size adjustment step includes: Performing a border removal operation on each of the 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 a 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; Wherein, the projection slice in the printing light intensity adjustment step is the scaled projection slice.

4. A volume bioprinting control method according to claim 3, characterized in that: The image scaling algorithm is a bi-triple interpolation algorithm.

5. A volume bioprinting control method according to claim 1, characterized in that: In the step S4, 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. A volume bioprinting control method according to claim 1, characterized in that: The step S32 specifically includes: S321, 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; S322, performing point-to-point correspondence multiplication of the discretized ramp filter and the undetermined window function with the same number of discretized sampling points to perform a windowing operation; S323, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head; S324, 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 excess 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.

7. A volume bioprinting control method according to claim 1, characterized in that: In step S3206, 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) Among them, in S3207, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in S3208, the window function corresponding to the printing evaluation index closest to 1 among the printing evaluation indexes with different parameters is used as the optimal window function.

8. A volume bioprinting control method 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