Projection filtering method, device, apparatus and storage medium suitable for volumetric bioprinting projection

Abstract

The present invention proposes a projection filtering method suitable for volume bioprinting projection, including: slicing the three-dimensional model equidistantly along the Z axis; performing a 360-degree Laden transform on each of the transverse slices to obtain W lines of projection data for each transverse slice; performing a fast Fourier transform on each line of projection data of the projection slices at different angles; 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 after the fast Fourier transform, and finally performing a fast Fourier inverse transform on each line of projection data after filtering each of the projection slices, thereby obtaining a filtered projection slice. By implementing a projection filtering method suitable for volume bioprinting projection provided by an embodiment of the present invention, the influence of star-shaped artifacts in volume bioprinting is reduced, and the printing effect is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of 3D printing, and in particular relates to a projection filtering method, device, equipment and storage medium suitable for volumetric bioprinting projection. Background Art

[0002] Volumetric bioprinting is an application of 3D printing technology in the biological field. It uses biological materials or bioactive substances to create biological tissues or organs with specific structures and functions by stacking or jetting layer by layer.

[0003] Volumetric bioprinting technology helps to solve some challenges in the medical field, such as tissue defect repair, organ transplantation, etc. By using biomaterials such as bio-ink, cells, proteins and other biomolecules can be accurately printed into three-dimensional structures to reconstruct damaged or missing tissues.

[0004] 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".

[0005] In volume bioprinting projection, very strong star-shaped artifacts are generated during projection printing, which seriously affect the printing effect and solidify the original non-solidified area. In existing volume bioprinting projection, it is generally necessary to perform filtering operations on the projection slices (corresponding to transverse slices in the prior art), but the filtering operation used in the existing volume bioprinting projection does not affect the filtering effect, and star-shaped artifacts still exist, thereby affecting the printing effect. Summary of the invention

[0006] The purpose of the present invention is to propose a projection filtering method, device, equipment and storage medium suitable for volumetric bioprinting projection, which can effectively solve the above-mentioned technical problems existing in the prior art.

[0007] In order to achieve the above object, an embodiment of the present invention provides a projection filtering method suitable for volumetric bioprinting projection, the method comprising the steps of:

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

[0009] 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 (w-1)*θ of the transverse slice obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form projection slices of the same angle;

[0010] S3, 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;

[0011] S4, 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;

[0012] S5, 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;

[0013] Wherein, the step S4 specifically includes:

[0014] S41, performing discrete sampling on 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;

[0015] S42, 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;

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

[0017] S44, 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.

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

[0019] S101, 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;

[0020] S102, 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°;

[0021] S103, 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;

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

[0023] S105, 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;

[0024] S106, 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:

[0025]

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

[0027] S107, 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;

[0028] S108, 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 S103, and comparing the printing evaluation indicators of the different window functions obtained after steps S104 to S107, 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.

[0029] As an improvement to the above, in step S106, 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:

[0030]

[0031] Wherein, in S107, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in S108, 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.

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

[0033] As an improvement of the above, the height of the three-dimensional model satisfies: [H / h]=M+1; the highest point of the three-dimensional model is located in the topmost transverse slice.

[0034] As an improvement of the above solution, the window function refers to the Kaiser window function under different parameters β, and the Kaiser window function formula is shown in formula (3):

[0035]

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

[0037] As an improvement of the above scheme, the M1 transverse slice samples correspond to transverse slices at heights of the 50th, 100th, 150th, 100th, 150th, 200th, 250th, 300th, 350th, 400th and 450th layers respectively.

[0038] Another embodiment of the present invention corresponds to a printing effect evaluation device suitable for volumetric bioprinting projection, comprising:

[0039] A transverse slice generation module is used to perform equidistant transverse slices on the three-dimensional model along the Z axis to obtain M transverse slices; wherein each transverse slice obtained is a transverse slice binary image, M≥1;

[0040] The projection data conversion module is used to perform 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°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form a projection slice of the same angle;

[0041] A fast Fourier transform module, used for performing 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;

[0042] A filtering module, used to perform a windowing operation on the ramp filter using a window function, and use the ramp filter after the windowing operation to filter each line of projection data of each projection slice after fast Fourier transformation;

[0043] The inverse fast Fourier transform module is used to perform inverse fast Fourier transform on each line of projection data after filtering of each projection slice.

[0044] Inverse fast Fourier transform to obtain filtered projection slices;

[0045] Wherein, the filtering module specifically includes:

[0046] A discrete sampling unit, used for discretizing 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;

[0047] A windowing operation unit, used for performing point-by-point corresponding multiplication of the discretized ramp filter and the undetermined window function with the same number of discretized sampling points to perform a windowing operation;

[0048] A data translation operation unit, used for performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head;

[0049] A filtering operation unit is used to perform a filtering operation by point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after each projection slice is subjected to a fast Fourier transform, and to discard the excess points between the ramp filter after the windowing and translation operation and the points of each line of projection data after each projection slice is subjected to a fast Fourier transform.

[0050] Yet another embodiment of the present invention provides an electronic device, comprising a processor and a memory, wherein the memory is used to store a computer program, the computer program comprises program instructions, and the processor is configured to call the program instructions to execute the projection filtering method applicable to volumetric bioprinting projection as described in any of the above embodiments.

[0051] Yet another embodiment of the present invention provides a computer-readable storage medium, characterized in that the computer-readable storage medium stores a computer program, wherein the computer program includes program instructions, and when the program instructions are executed by a processor, the processor executes the projection filtering method applicable to volumetric bioprinting projection as described in any of the above embodiments.

[0052] Compared with the prior art, the embodiments of the present invention provide a projection filtering method, device, equipment and storage medium suitable for volumetric bioprinting projection. By adding Fourier transform, filtering and inverse Fourier transform processes in the filtered back-projection algorithm in volumetric bioprinting, the three-dimensional model is specifically sliced ​​equidistantly along the Z axis to obtain M transverse slices; a 360-degree Laden transform is performed on each of the transverse slices to obtain W lines of projection data for each transverse slice; a fast Fourier transform is performed on each line of projection data of the projection slices at different angles of 360 degrees on the side of the three-dimensional model; a window function is used to perform a windowing operation on the ramp filter, and each line of projection data of each projection slice after the fast Fourier transform is filtered using the ramp filter after the windowing operation, including filtering the ramp filter. The slope filter performs discretization sampling in the frequency domain, performs point-to-point multiplication of the discretized slope filter and the undetermined window function with the same number of sampling points discretized to perform a windowing operation, performs a translation operation on the windowed slope filter data, thereby moving the last half of the slope filter data to the head, performs point-to-point multiplication of the slope filter after windowing and translation operation and each line of projection data after each projection slice is subjected to fast Fourier transformation to perform a filtering operation, and discards the extra points of the slope filter after windowing and translation operation compared with each line of projection data after each projection slice is subjected to fast Fourier transformation; finally, each line of projection data of each projection slice is subjected to fast Fourier inverse transformation after filtering, thereby obtaining a filtered projection slice. After practical verification, by implementing a projection filtering method suitable for volume bioprinting projection provided by an embodiment of the present invention, the influence of star-shaped artifacts in volume bioprinting is reduced, and the printing effect is improved. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0054] Figure 1 A schematic flow chart of a projection filtering method suitable for volumetric bioprinting projection provided in an embodiment of the present invention.

[0055] Figure 2 A schematic structural diagram of a transverse slice generated by a projection filtering method suitable for volumetric bioprinting projection provided by an embodiment of the present invention.

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

[0057] Figure 4 The process of generating (transforming) the obtained projection slices in the printing method evaluated by a projection filtering method suitable for volumetric bio-printing projection provided by an embodiment of the present invention is demonstrated.

[0058] Figure 5 A schematic diagram of volumetric printing of a three-dimensional model of the obtained projection slices in a printing method evaluated by a projection filtering method suitable for volumetric bio-printing projection provided by an embodiment of the present invention.

[0059] Figure 6 It is a schematic diagram of the unfiltered projection slice of the ring wave model at 0 degrees.

[0060] Figure 7 It is a schematic diagram of the filtered projection slice of the ring wave model after executing steps S3 to S5 at 0 degrees.

[0061] Figure 8 The figure shows a comparison between a binary image of a transverse slice corresponding to the 400th transverse slice of the annular wave 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).

[0062] 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 β.

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

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

[0065] Fig.12 A structural block diagram of a printing effect evaluation device suitable for volumetric bioprinting projection provided in an embodiment of the present invention.

[0066] Fig.13 It is a structural schematic diagram of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0068] Figure 1 A flowchart of a projection slicing method suitable for volume bioprinting provided in an embodiment of the present application, the method comprising steps S1 to S5:

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

[0070] 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 (w-1)*θ of the transverse slice obtained by Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form projection slices of the same angle;

[0071] S3, 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;

[0072] S4, 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;

[0073] S5, 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;

[0074] Wherein, the step S4 specifically includes:

[0075] S41, performing discrete sampling on 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;

[0076] S42, 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;

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

[0078] S44, 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.

[0079] Below, each step of a projection filtering method suitable for volumetric bioprinting projection provided in an embodiment of the present application will be described in detail.

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

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

[0082] Specifically, the following is combined Figure 2 , how to obtain the horizontal slice 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 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:

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

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

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

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

[0087]

[0088]

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

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

[0091]

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

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

[0094] 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 450th transverse slice of the three-dimensional model are respectively shown.

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

[0096] 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.01mm.

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

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

[0099] It can be understood that after the transverse slice binary images of all transverse slices of the three-dimensional model are obtained through step S1, the embodiment of the present invention does not directly use the transverse slices as projection slices for 3D volume printing, but performs a 360-degree Radon transform on all transverse slices of the three-dimensional model to obtain 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 printing in angular order.

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

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

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

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

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

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

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

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

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

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

[0110] Specifically, in the above step S2, 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 Projected slices at 0° are shown.

[0111] like Figure 5As shown, the projection slices obtained by the projection filtering method suitable for volumetric bio-printing provided by the present invention are used to perform volumetric printing of a three-dimensional model. During the volumetric 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.

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

[0113] Returning to the projection filtering method applicable to volume bioprinting projection in the embodiment of the present invention, Figure 6 to Figure 8 Steps S3 to S5 are described in detail. Figure 6 It is a schematic diagram of the unfiltered projection slice of the ring wave model at 0 degrees. Figure 7 It is a schematic diagram of the filtered projection slice of the ring wave model after executing steps S3 to S5 at 0 degrees. Figure 8 The figure shows a comparison between a binary image of a transverse slice corresponding to the 400th transverse slice of the annular wave 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).

[0114] Specifically, the projection filtering method for volume bioprinting projection provided by the embodiment of the present invention is to innovatively add 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 star 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 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, 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 artifacts, and improve image clarity.

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

[0116] By comparison Figure 6 and Figure 7 It can be seen that Figure 6 What is shown is a 90-degree projection slice of the ring wave model obtained by executing the above steps S1 and S2 (i.e., without executing steps S3 to S5), with a height of 500 pixels and a width of 800 pixels. Figure 7 Will Figure 6 The 0-degree unfiltered projection slice of the ring wave model shown is a 0-degree filtered projection slice of the ring wave model after the filtering of steps S3 to S5 is performed by the Kaiser window function with β=6, Figure 7 It can be seen that the edges are clear and the star-shaped artifacts are effectively reduced.

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

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

[0119]

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

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

[0122]

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

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

[0125] S101, 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;

[0126] S102, 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°;

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

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

[0129] S105, 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;

[0130] S106, 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:

[0131]

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

[0133] S107, 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;

[0134] S108, 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 S103, and comparing the printing evaluation indicators of the different window functions obtained after steps S104 to S107, 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.

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

[0136] It can be understood that in this embodiment, the M1 transverse slice samples used as sampling are the ring wave 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.

[0137] It can be understood that in step S2, each angle of the 360 ​​degrees of the 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.

[0138] In this embodiment, after the projection data of each angle of each transverse slice sample is obtained in step S102, a filter back-projection operation is performed on the projection data of each angle of each transverse slice in step S103, 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.

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

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

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

[0142] S1033, 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.

[0143] Wherein, the step S32 specifically includes:

[0144] S10321, 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;

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

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

[0147] S10324, 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.

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

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

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

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

[0152]

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

[0154] Furthermore, in step S104 to step S105, 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.

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

[0156] Next, combine Figures 9 to 11 , the steps S106 to S108 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 represents the printing evaluation index. Specifically, Fig.11 The figure shows the change of the printing evaluation index (average evaluation value) of the Kaiser window function with β of 4 with the critical curing light intensity value.

[0157] Specifically, in steps S106 to S107, 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.

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

[0159]

[0160] Wherein, when executing the solution of the above formula (2), in the step S107, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in the step S108, 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.

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

[0162] Then, in step S108, 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 S104 to S107, 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 best window function.

[0163] like Fig.12 As shown, an embodiment of the present invention provides a printing effect evaluation device suitable for volumetric bioprinting projection, which includes:

[0164] The transverse slice generation module 110 is used to perform equidistant transverse slices on the three-dimensional model along the Z axis to obtain M transverse slices; wherein each transverse slice obtained is a transverse slice binary image, and M≥1;

[0165] The projection data conversion module 120 is used to perform 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°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form a projection slice of the same angle;

[0166] A fast Fourier transform module 130 is used to perform 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;

[0167] A filtering module 140 is used to perform a windowing operation on the shelving filter using a window function, and use the shelving filter after the windowing operation to filter each line of projection data of each projection slice after fast Fourier transformation;

[0168] A fast Fourier inverse transform module 150 is used to perform a fast Fourier inverse transform on each line of projection data of each projection slice after filtering, so as to obtain a filtered projection slice;

[0169] Wherein, the filtering module specifically includes:

[0170] A discrete sampling unit, used for discretizing 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;

[0171] A windowing operation unit, used for performing point-by-point corresponding multiplication of the discretized ramp filter and the undetermined window function with the same number of discretized sampling points to perform a windowing operation;

[0172] A data translation operation unit, used for performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head;

[0173] A filtering operation unit is used to perform a filtering operation by point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after each projection slice is subjected to a fast Fourier transform, and to discard the excess points between the ramp filter after the windowing and translation operation and the points of each line of projection data after each projection slice is subjected to a fast Fourier transform.

[0174] The specific implementation of the printing effect evaluation device applicable to volumetric bioprinting projection in this embodiment can refer to the description of the projection filtering method applicable to volumetric bioprinting projection in the above-mentioned embodiment, which will not be repeated here.

[0175] like Fig.13 As shown, an embodiment of the present invention provides an electronic device 300, including a memory 310 and a processor 320, wherein the memory 310 is used to store one or more computer instructions, and the processor 320 is used to call and execute the one or more computer instructions, so as to implement any of the above-mentioned projection filtering methods applicable to volumetric bioprinting projection.

[0176] That is, the electronic device 300 includes: a processor 320 and a memory 310, in which computer program instructions are stored, wherein when the computer program instructions are executed by the processor, the processor 320 executes any of the above-mentioned projection filtering methods applicable to volumetric bioprinting projection.

[0177] Furthermore, if Fig.13 As shown, the electronic device 300 further includes a network interface 330 , an input device 340 , a hard disk 350 , and a display device 360 ​​.

[0178] The above-mentioned interfaces and devices can be interconnected through a bus architecture. The bus architecture can be a bus and bridge that can include any number of interconnected buses. Specifically, one or more central processing units (CPUs) represented by processor 320 and various circuits of one or more memories represented by memory 310 are connected together. The bus architecture can also connect various other circuits such as peripheral devices, voltage regulators, and power management circuits together. It can be understood that the bus architecture is used to achieve connection and communication between these components. In addition to the data bus, the bus architecture also includes a power bus, a control bus, and a status signal bus, which are all well known in the art, so they will not be described in detail herein.

[0179] The network interface 330 can be connected to a network (such as the Internet, a local area network, etc.), obtain relevant data from the network, and save it in the hard disk 350.

[0180] The input device 340 can receive various instructions input by the operator and send them to the processor 320 for execution. The input device 340 can include a keyboard or a pointing device (for example, a mouse, a trackball, a touch pad or a touch screen, etc.).

[0181] The display device 360 ​​can display the result obtained by the processor 320 executing the instruction.

[0182] The memory 310 is used to store programs and data necessary for the operation of the operating system, as well as data such as intermediate results during the calculation process of the processor 320.

[0183] It is understood that the memory 310 in the embodiment of the present invention may be a volatile memory or a non-volatile memory, or may include both volatile and non-volatile memories. Among them, the non-volatile memory may be a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or a flash memory. The volatile memory may be a random access memory (RAM), which is used as an external cache. The memory 310 of the apparatus and method described herein is intended to include, but is not limited to, these and any other suitable types of memory.

[0184] In some implementations, the memory 310 stores the following elements, executable modules or data structures, or a subset thereof, or an extended set thereof: an operating system 311 and application programs 312 .

[0185] The operating system 311 includes various system programs, such as a framework layer, a core library layer, a driver layer, etc., which are used to implement various basic services and process hardware-based tasks. The application 312 includes various application programs, such as a browser, etc., which are used to implement various application services. The program for implementing the method of the embodiment of the present invention can be included in the application 312.

[0186] The processor 320, when calling and executing the application and data stored in the memory 310, specifically, the program or instruction stored in the application 312, slices the three-dimensional model horizontally along the Z axis at equal intervals to obtain M horizontal slices; performs a 360-degree Laden transform on each of the horizontal slices to obtain W lines of projection data for each horizontal slice; performs 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; performs a windowing operation on the ramp filter using a window function, and uses the ramp filter after the windowing operation to filter each line of projection data of each of the projection slices after the fast Fourier transform, including discretizing the ramp filter in the frequency domain. The discretized ramp filter is point-to-point multiplied with the undetermined window function discretized with the same number of sampling points to perform a windowing operation, the windowed ramp filter data is translated to move the last half of the ramp filter data to the head, the ramp filter after the windowing and translation operation is point-to-point multiplied with each line of projection data of each projection slice after fast Fourier transform to perform a filtering operation, and the extra points of the ramp filter after the windowing and translation operation compared with the points of each line of projection data of each projection slice after fast Fourier transform are discarded; finally, each line of projection data of each projection slice after filtering is inversely transformed by fast Fourier transform to obtain a filtered projection slice.

[0187] The projection filtering method for volume bioprinting projection disclosed in the above embodiment of the present invention can be applied to the processor 320, or implemented by the processor 320. The processor 320 may be an integrated circuit chip with signal processing capabilities. In the implementation process, each step of the above method can be completed by the hardware integrated logic circuit or software instructions in the processor 320. The above processor 320 can be a general processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, and can implement or execute the methods, steps and logic block diagrams disclosed in the embodiments of the present invention. The general processor can be a microprocessor or the processor can also be any conventional processor, etc. The steps of the method disclosed in conjunction with the embodiment of the present invention can be directly embodied as a hardware decoding processor to execute, or can be executed by a combination of hardware and software modules in the decoding processor. The software module can be located in a mature storage medium in the field such as a random access memory, a flash memory, a read-only memory, a programmable read-only memory or an electrically erasable programmable memory, a register, etc. The storage medium is located in the memory 310, and the processor 320 reads the information in the memory 310 and completes the steps of the above method in combination with its hardware.

[0188] It is understood that the embodiments described herein can be implemented in hardware, software, firmware, middleware, microcode or a combination thereof. For hardware implementation, the processing unit can be implemented in one or more application specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), general purpose processors, controllers, microcontrollers, microprocessors, other electronic units for performing the functions described in the present application or a combination thereof.

[0189] For software implementation, the techniques described herein can be implemented by modules (e.g., procedures, functions, etc.) that perform the functions described herein. The software code can be stored in a memory and executed by a processor. The memory can be implemented in the processor or outside the processor.

[0190] Specifically, the processor 320 is also used to read the computer program and execute any of the above-mentioned projection filtering methods applicable to volumetric bioprinting projection.

[0191] The present application also provides a computer-readable storage medium, which stores a computer program. The computer program includes program instructions. When the program instructions are executed by a processor, the processor executes the above method, such as executing the method executed by the above electronic device, which is not repeated here.

[0192] Optionally, the storage medium involved in the present application, such as computer-readable storage medium, may be non-volatile or volatile.

[0193] Optionally, the computer-readable storage medium may mainly include a program storage area and a data storage area, wherein the program storage area may store an operating system, an application required for at least one function, etc.; the data storage area may store data created according to the use of blockchain nodes, etc. Among them, the blockchain referred to in this application is a new application mode of computer technologies such as distributed data storage, point-to-point transmission, consensus mechanism, encryption algorithm, etc. Blockchain is essentially a decentralized database, a string of data blocks generated by cryptographic methods, each of which contains a batch of network transaction information, which is used to verify the validity of its information (anti-counterfeiting) and generate the next block. The blockchain may include the underlying blockchain platform, the platform product service layer, and the application service layer.

[0194] It should be noted that, for the above-mentioned various method embodiments, for the sake of simplicity of description, they are all expressed as a series of action combinations, but those skilled in the art should be aware that this application is not limited by the order of the actions described, because according to this application, some steps can be performed in other orders or simultaneously. Secondly, those skilled in the art should also be aware that the embodiments described in the specification are all preferred embodiments, and the actions and modules involved are not necessarily required by this application.

[0195] In the several embodiments provided in the present application, it should be understood that the disclosed methods and devices can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0196] In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may be physically included separately, or two or more units may be integrated into one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of hardware plus software functional units.

[0197] The above-mentioned integrated unit implemented in the form of a software functional unit can be stored in a computer-readable storage medium. The above-mentioned software functional unit is stored in a storage medium, including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform some steps of the sending and receiving method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (Read-Only Memory, referred to as ROM), random access memory (Random Access Memory, referred to as RAM), disk or optical disk and other media that can store program codes.

[0198] 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 projection filtering method suitable for volumetric bioprinting projection, characterized in that: The method comprises the steps of Steps: 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°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form projection slices of the same angle; S3, 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; S4, use the window function to perform a windowing operation on the ramp filter, and use the ramp filter after the windowing operation to filter each image. The projection slice is filtered for each line of projection data after fast Fourier transformation; S5, performing inverse fast Fourier transform on each line of projection data after filtering of each projection slice, so as to Get the filtered projection slice; Wherein, the step S4 specifically includes: S41, performing discrete sampling on 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; S42, performing point-by-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; S43, performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head; S44, 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.

2. A projection filtering method suitable for volumetric bioprinting projection according to claim 1, characterized in that In, The window function is an optimal window function, and the optimal window function is determined by the following steps: S101, 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; S102, 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 the transverse slice sample at an angle of (w-1)*θ of the Laden transformation, w=1, 2...W, 0°<θ≤1°, W*θ=360°; S103, 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; S104, normalizing the pixel value of each back-projection reconstructed image to an integer between 0 and 255, to obtain M1 normalized reconstructed images; S105, 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; S106, 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; is the value of the pixel at the position (x, y) with the upper left corner as the coordinate (1, 1) and the positive direction downward and left in the horizontal slice sample; is the value of the pixel at the position (x, y) with the upper left corner as the coordinate (1, 1) and the positive direction downward and left in the simulated printed horizontal slice; S107, 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; S108, 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 S103, and comparing the printing evaluation indicators of the different window functions obtained after steps S104 to S107, 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.

3. A projection filtering method suitable for volumetric bioprinting projection according to claim 2, characterized in that: In step S106, 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 S107, the average value closest to 1 is used as the printing evaluation index of the undetermined window function; in S108, 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.

4. A projection filtering method suitable for volumetric bioprinting projection 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.

5. A projection filtering method suitable for volumetric bioprinting projection according to claim 1, characterized in that: The height of the three-dimensional model satisfies: ; The highest point of the three-dimensional model is located in the topmost transverse slice.

6. A projection filtering method suitable for volumetric bioprinting projection according to claim 1, characterized in that: The window function refers to the Kaiser window function under different parameters β. The Kaiser window function formula is shown in formula (3): Formula (3) Where 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.

7. A projection filtering method suitable for volumetric bioprinting projection according to claim 2, characterized in that: The transverse slice samples in M1 correspond to transverse slices at heights of the 50th, 100th, 150th, 100th, 150th, 200th, 250th, 300th, 350th, 400th and 450th layers respectively.

8. A printing effect evaluation device suitable for volumetric bioprinting projection, characterized in that: include: A transverse slice generation module is used to perform equidistant transverse slices on the three-dimensional model along the Z axis to obtain M transverse slices; wherein each transverse slice obtained is a transverse slice binary image, M≥1; The projection data conversion module is used to perform 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°; the projection data of the same angle of the M transverse slices are stacked in sequence along the direction of the Z axis to form projection slices of the same angle; A fast Fourier transform module, used for performing 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; The filtering module is used to perform a windowing operation on the ramp filter using a window function, and to use the ramp filter after the windowing operation The device filters each line of projection data of each projection slice after fast Fourier transformation; The inverse fast Fourier transform module is used to perform inverse fast Fourier transform on each line of projection data after filtering of each projection slice. Inverse fast Fourier transform to obtain filtered projection slices; Wherein, the filtering module specifically includes: A discrete sampling unit, used for discretizing 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; A windowing operation unit, used for performing point-by-point corresponding multiplication of the discretized ramp filter and the undetermined window function with the same number of discretized sampling points to perform a windowing operation; A data translation operation unit, used for performing a translation operation on the windowed shelving filter data, thereby moving the last half of the shelving filter data to the head; A filtering operation unit is used to perform a filtering operation by point-to-point multiplication of the ramp filter after the windowing and translation operation with each line of projection data after each projection slice is subjected to a fast Fourier transform, and to discard the excess points between the ramp filter after the windowing and translation operation and the points of each line of projection data after each projection slice is subjected to a fast Fourier transform.

9. An electronic device, characterized in that: It includes a processor and a memory, wherein the memory is used to store a computer program, the computer program includes program instructions, and the processor is configured to call the program instructions to execute the projection filtering method suitable for volumetric bioprinting projection as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, which includes program instructions. When the program instructions are executed by a processor, the processor executes the projection filtering method suitable for volumetric bioprinting projection as described in any one of claims 1 to 7.

Images