JP7197316B2 - Image processing method, computer program and recording medium - Google Patents

Description

The present invention relates to an image processing technique for determining the presence or absence of pulsation from an image of cultured myocardial cell colonies.

In fields such as regenerative medicine and drug discovery, for example, iPS cells (induced Pluripotent Stem Cells) and ES cells (Embryonic
Techniques for artificially producing cardiomyocytes by differentiating and culturing pluripotent stem cells such as stem cells are being studied. Since well-cultured cardiomyocytes exhibit periodic beating motion, the presence or absence of beating can be determined by observing samples containing cardiomyocyte colonies under a microscope or the like in order to evaluate the culture state. done. Since the results of visual judgment vary depending on the subjectivity of the observer, and observing a large number of colonies is a burden on the observer, the presence or absence of pulsation is automatically detected by data analysis processing of the image of the sample. A method to do so is proposed.

For example, in the technique described in Patent Document 1, a motion vector is detected from time-series images captured at predetermined time intervals, and among the regions in which the motion vector is detected, those in which a significant change in the internal structure is observed are pulsations. Detected as a region.

JP 2014-075999 A

In a sample at the initial stage of differentiation culture, a large number of small colonies formed by cardiomyocytes differentiated from stem cells are dispersed in the sample. More specifically, some cardiomyocytes aggregate spherically to form cell colonies, also called spheres, and such colonies are scattered within the sample. Those colonies beat in mutually independent cycles, and colonies that do not beat may also be included. In addition, the movement of the colony during pulsation is not necessarily contraction/relaxation, and often appears as periodic vibration in the image. Furthermore, in images of specimens in which colonies are present at a high density, the colonies are often in contact with each other.

The prior art described above is suitable for detecting rhythmic beats in relatively large colonies. However, as described above, in a sample in which a large number of small colonies are scattered and individual colonies beat individually, it is not necessarily suitable for determining the presence or absence of individual colony beats from the captured image. Conceivable. Therefore, it is desired to establish an image processing technique for pulsation detection that is suitable when there are many such small colonies.

The present invention has been made in view of the above problems, and an object of the present invention is to provide an image processing technique suitable for determining the presence or absence of pulsation in individual colonies from an image containing a plurality of colonies of cultured cardiomyocytes. and

In order to achieve the above object, an image processing method according to the present invention provides a step of acquiring a plurality of original images obtained by imaging a sample containing colonies of cultured cardiomyocytes at different imaging times; a step of extracting a colony region corresponding to the colony from each of the first image and the second image; detecting a displacement of the colony region between the first image and the second image; and based on the detection result. and determining whether or not the colony pulsates.

In one aspect of the present invention, for one target colony region among the colony regions, the magnitude of the difference between the first image and the second image within a predetermined search range is the displacement of the target colony region. The search range is an annular area obtained by expanding the outline of the target colony area in the first image to a predetermined width, and the outline of the target colony area in the second image to the predetermined width. This area corresponds to the logical sum with the expanded annular area.

In another aspect of the present invention, for one target colony region among the colony regions, the magnitude of the difference between the first image and the second image within a predetermined search range is calculated as the target colony region. and the search range corresponds to the exclusive OR of the internal area of the target colony area in the first image and the internal area of the target colony area in the second image. It is an annular region that includes at least the region.

In the invention configured as described above, the colony area is extracted in each of the first image and the second image captured at different times, and the displacement of the colony area between the two images is detected. This displacement is detected as the difference between the colony regions extracted in the first image and the second image. However, instead of a simple difference between colony regions, a difference within a search range set based on the contours of the extracted colony regions is obtained.

Although details will be described later, in both of the above-described two modes, a search range is set around the contours of the colony areas detected in the first image and the second image. Therefore, the magnitude of the calculated difference strongly reflects the magnitude of movement of the contour portion of the colony region. On the other hand, the displacement inside or outside the colony, which is far from the outline portion and is not included in the search range, does not affect the value of the difference. Therefore, the size of the difference obtained for one search range represents the size of displacement of one colony.

According to the findings of the inventors of the present application, for small colonies such as cardiomyocyte spheres, the displacement of each colony in the sample can be estimated by using the magnitude of the difference obtained in this way as an index of the magnitude of displacement of the colony. It can be grasped individually and accurately. Then, it is possible to determine the presence or absence of pulsation of the colony from the magnitude of the displacement grasped in this way. Even if there are many colonies in the original image, only the displacement around the contour of each colony is detected, so the presence or absence of pulsation can be determined without being affected by the behavior of other colonies. It is possible.

The image processing method according to the present invention can be executed by a computer. In this sense, the present invention can be implemented as a computer program for causing a computer to execute the above processes. Moreover, it is also possible to implement it as a recording medium recording the computer program.

As described above, according to the present invention, the image difference around the contour of each colony identified from an image containing a plurality of colonies of cultured myocardial cells is used as an index of the displacement of the colony. Therefore, even if an image contains a large number of colonies, the presence or absence of pulsation of individual colonies can be appropriately determined without being affected by the behavior of other colonies.

4 is a flowchart showing image processing according to the embodiment; FIG. 4 is a diagram schematically illustrating one frame image; It is a flow chart which shows a pretreatment process. It is a flow chart which shows a sphere extraction process. FIG. 10 is a diagram for explaining an image edge correction process; It is a figure which shows the example of the image after edge detection processing. FIG. 4 is a diagram showing an example of a contour area image; It is a figure which shows an example of a binarized image. It is a figure which shows an example of a logical product composite image. It is a figure which shows an example of a distance image. FIG. 10 is a diagram showing the concept of unnecessary object elimination based on distance transform results; FIG. 4 is a diagram showing an example of a central image; FIG. 4 is a diagram showing the concept of object division; FIG. 4 is a diagram showing the correspondence of spheres between different frames; FIG. 4 is a diagram showing the principle of thinning processing; 4 is a flow chart showing a motion detection process; FIG. 4 is a diagram showing the relationship between a reference image and a comparison image; It is an example of an image showing a difference between a reference image and a comparative image. It is a figure which shows the 1st aspect of the search range setting method. It is a figure which shows the 2nd aspect of the search range setting method. FIG. 10 is a diagram showing a logical product synthesized image of an eroded and expanded image and a search range image; It is a flow chart which shows a pulsation judging process. It is a figure which shows the determination example of a motion direction. It is a configuration example of a computer device that executes the image processing of this embodiment.

<Description of overview>
An embodiment of the image processing method according to the present invention will be described below. The image processing in this embodiment is processing for detecting cell colonies from an image of a sample containing cultured myocardial cell colonies and determining whether or not each cell colony is beating. For example, iPS cells (induced Pluripotent Stem Cells) and ES cells (Embryonic
A sample containing colonies of myocardial cells differentiated and cultured from pluripotent stem cells such as Stem Cells can be used as an object to be processed in this embodiment.

At the initial stage of such culture, several cardiomyocytes differentiated from stem cells aggregate to form small colonies called roughly spherical spheres. A sample can contain a large number of spheres that beat independently of each other and some do not. In addition, the sample may also contain undifferentiated cells, waste products excreted in the course of cell metabolism, and the like. The image processing of this embodiment, which will be described below, aims to detect individual spheres from such an image and individually determine the presence or absence of pulsation.

This image processing is performed on images captured by an imaging device having a moving image capturing function that repeatedly captures images of a sample containing myocardial cell spheres at predetermined time intervals. The image processing may be performed as one of the post-processing of the image data after the imaging by the imaging device, and the computer device that receives the image data corresponding to the captured image from the imaging device or appropriate storage means processes the main image. processing may be performed.

Here, a mode in which a general-purpose computer device controls a known imaging device to acquire an image and executes the image processing of the present embodiment on the image will be described as an example. The configuration of the imaging device is not particularly limited, and any device may be used as long as it has a function of capturing a moving image of a sample containing cultured cells or the like together with a culture medium and outputting it as digital image data. For example, an imaging device in which a digital video camera is combined with a microscope can be used. A bright field image, for example, can be preferably used as the image. As for the computer device, for example, a computer having a general hardware configuration commercialized as a personal computer can be used.

FIG. 1 is a flowchart showing image processing according to this embodiment. First, the outline of each step of image processing in this embodiment will be described with reference to FIG. 1, and then the specific contents of the processing in each step will be described in detail. Note that this image processing includes two processes: a sphere detection process (steps S101 to S103) for detecting spheres from the captured image, and a pulsation detection process (steps S105 to S107) for detecting the pulsation of each detected sphere. It can be roughly divided into two types of processing. In addition, a thinning process is provided between these to achieve a frame rate suitable for each process.

The computer device captures a moving image of a sample containing colonies of myocardial cells cultured in a medium at an appropriate frame rate (step S101). As for the frame rate, it is preferable that the imaging cycle is sufficiently shorter than the cycle assumed to be the cardiac muscle cell beating cycle. can do. In addition, the duration of imaging can be set to a time, for example, several seconds, during which the movement of the pulsating cell colony can be reliably captured over one or more cycles.

FIG. 2 is a diagram schematically illustrating one frame image. In the frame image If, a large number of objects corresponding to spheres Sp, which are formed by gathering cardiomyocytes in a substantially spherical shape, have a higher density than the background in a substantially uniform and low-density (high-brightness) background BG. distributed. In addition, objects corresponding to objects De other than spheres, such as undifferentiated cells, free cells, and clumps of waste excreted from cells, also exist in the image. If the density relationship between the background and the object is opposite to that described above, the gradation in the following description should also be read in reverse. In order to make the drawing easier to see, not all objects are labeled, but only some typical objects are labeled. The same applies to each of the following figures.

Multiple objects may be in contact with each other in the image, as exemplified by reference numerals AD in FIG. The reason for such coupling is that the objects are actually touching each other in the sample, or they are separated in the sample but appear to be touching in the image.

A preprocessing step for correcting minute disturbances in the image is executed for each frame of the image captured in this way (step S102). Then, a sphere extraction step is performed on the preprocessed frame image (step S103). This processing is processing for extracting what appears to be a sphere from among the objects scattered in the image and specifying the contour position of the object. This identifies which area the sphere occupies in the image. If the image contains multiple spheres, their contours are identified individually.

Subsequently, frame thinning is performed (step S104). As will be described later, the pulsation detection of the sphere in this embodiment is performed based on the magnitude of displacement of the sphere between frame images. Here, imaging of the sample is performed at a sufficiently short period with respect to the pulsation period of the sphere, but the displacement of the sphere is small between images with a too small difference in imaging time, making it difficult to detect it. Therefore, it is not always realistic to try to detect the motion of the sphere using all frame images.

On the other hand, if the difference in imaging time is too large, there is a risk of overlooking the periodic displacement. By thinning out the images at appropriate intervals in consideration of the pulsation cycle, it becomes possible to more reliably capture the displacement of the sphere. It also contributes to greatly reducing the amount of calculation for image processing. For this purpose, frame thinning is performed to select a part of the picked-up frame images at a predetermined sampling interval.

Next, motion of the sphere between images is detected (step S105). As described above, the movement of the sphere is detected from the displacement of the sphere between two images captured at different times. Then, the presence or absence of pulsation of each sphere is determined from the amount of displacement thus detected (step S106). Displacement of an object in an image is not limited to pulsation. The presence or absence of pulsation is determined based on whether the detected displacement is caused by pulsation. The determination result is output in an appropriate format and presented to the user (step S107). The above is the overall flow of image processing according to this embodiment. Specific contents of each process will be described below in sequence.

<Sphere detection processing>
FIG. 3 is a flow chart showing the pretreatment process. Contrast correction (step S201) and noise removal (step S202) are performed on the captured image of each frame in order to make the conditions for image processing, which will be described later, constant. Known processes can be applied as these processes, and the description thereof is omitted here. For example, a spatial filtering technique such as a bilateral filter can be used as noise removal processing. Depending on the imaging conditions, these preprocessing may be omitted.

Next, the sphere extraction step (step S103) will be described. In the sphere extraction step, regions of various objects existing in the image are identified, spheres Sp are selected from among them, and the spheres Sp that are combined in the image are separated by appropriate boundaries. , to identify objects in the image that correspond to individual spheres. This processing is executed for all frame images If obtained by imaging.

FIG. 4 is a flow chart showing the sphere extraction process. First, one frame image is selected (step S301), and image edge correction is performed on the selected frame image If (step S302). As shown in FIG. 2, an object in the image near the edge of the frame image If may be partly over the edge of the image and not entirely within the image. In the processing to be described later, there is a step of extracting an area surrounded by a closed curve that can be regarded as a contour, but in order to avoid missing an object whose contour is interrupted at the edge of the image, image edge correction is performed. .

FIG. 5 is a diagram for explaining the image edge correction process. As shown in FIG. 5A, a frame-shaped region Rf having a width of one pixel is provided so as to surround the peripheral portion of the frame image If. The pixel value (luminance value) of this area is set to be approximately the same as or lower in density (higher luminance) than the background BG. For example, the density can be 0 (ie, white with maximum brightness). This region Rf may be provided by converting the pixel values of pixels corresponding to the edges of the image to the above values, and adding virtual pixels having the above pixel values outside the edges of the image. may be provided by

In the processed frame image If, as shown in FIG. 5(b), all the pixels at the edges of the image have the same or lower density than the background density. For this reason, objects (for example, indicated by symbols A and B in the figure) that overlap the edges of the image also have closed contours in the image. This avoids the omission of object extraction due to the discontinuity of the contour. Note that in FIG. 5B, the peripheral edge of the frame image If is indicated by a dotted line in order to make it easier to see the edge of the image.

Returning to FIG. 4, edge detection processing is executed to detect portions where luminance values change significantly on the frame image If whose image edges have been corrected in this way (step S303). As the edge detection processing, a known method such as spatial filtering represented by a Sobel filter can be applied.

FIG. 6 is a diagram showing an example of an image after edge detection processing. In the image (edge image) Ie after edge detection processing corresponding to the frame image If shown in FIG. 2, the boundary between the object and the background in the image, that is, the portion corresponding to the outer periphery of the object is emphasized and detected as an edge. . Edges due to internal textures of the object may also be included, as exemplified by the letter A in the figure. In addition, since the image edge correction is performed in advance, as exemplified by the symbol B in the figure, a closed contour is also formed for the object at the image edge. Appropriate correction processing, such as expansion and contraction processing, may be further performed to compensate for minute omissions in the detected contour.

When the edges thus detected form a closed curve, the closed curve is highly likely to be the contour of an object in the image. Therefore, by extracting the inner region of the contour surrounded by the closed curve in the edge image Ie and painting it (step S304), the inside and the outside of the contour are distinguished. Here, the region inside the filled contour is called "contour region", and the image after filling is called "contour region image".

FIG. 7 is a diagram showing an example of a contour area image. This contour area image Ir corresponds to an image in which the area surrounded by the closed curve in the edge image Ie of FIG. Here, the colors are divided according to the example of general image processing, but the shading may be reversed. The same applies to the following processing.

In the contour area image Ir, the inside of the object surrounded by the closed curve is completely painted out, so the information corresponding to the internal texture of the object is lost. Also, when compared with the original frame image If shown in FIG. 2, as illustrated with reference numerals A and B in FIG. Some areas will be filled in. That is, the background area, which should be outside the object, is treated as inside the object.

Actions are then taken to avoid this problem. That is, an image is generated by binarizing the original frame image If using an appropriate threshold value (step S305), and a composite image is generated by logically synthesizing the binarized image and the contour area image Ir (step S305). S306). As the frame image If to be binarized, one after preprocessing is preferably used. Note that the threshold value for binarization may not be a fixed value, but may be changed as a user setting parameter, for example.

FIG. 8 is a diagram showing an example of a binarized image. More specifically, in FIG. 8, the image obtained by preprocessing the frame image If in FIG. is shown in a dark color (black). In the binarized image Ib, the contour portion of the object and the high-density area inside it are represented in white, while the background area and some low-density areas inside the object, which were low-density in the original frame image If, are represented in white. is represented in black. As illustrated with symbols A and B in FIG. 8, the background portion surrounded by objects is also represented in black. Note that the white areas appearing in the upper left corner and upper right corner of image Ib, which are originally the background, symbolically indicate that background density unevenness may affect the result of binarization. It is a thing.

FIG. 9 is a diagram showing an example of a logical product synthesized image. If the contour area image Ir shown in FIG. 7 and the binarized image Ib shown in FIG. , and other areas are represented in black. As a result, the background area surrounded by the object, which was filled in white in the outline area image Ir of FIG. 7, is represented in the same black color as the other background portions. On the other hand, the black area due to the internal texture of the object contained in the binarized image Ib is also preserved in the synthesized image Is.

Therefore, the synthesized image Is shown in FIG. 9 is more accurate than the outline region image Ir shown in FIG. It can be said that it is effective. In addition, the synthesized image Is is more accurate than the binarized image Ib in FIG. 8 in the sense that the contours detected as clear edges define the boundary between the object and the background without being affected by the density unevenness of the background. It can be said that it is effective. In other words, the composite image Is combines the advantages of the contour extraction-based method and the binarization-based method for the purpose of separating the object and the background in the image.

Through the processing up to this point, the area occupied by the object in the image and the background area are accurately distinguished. However, as described above, the objects in the image include unnecessary objects other than the sphere, which is the object of processing in this embodiment, and there are cases where a plurality of objects are in contact with each other and combined in the image. Therefore, processing for removing unnecessary objects and dividing combined objects is required.

For this purpose, a structural analysis of each object in the image is performed. Specifically, distance conversion is performed on the synthesized image Is generated in step S306 (step S307). Distance transform in a two-dimensional binary image is general image processing that can be executed even with general-purpose image processing software, and represents each pixel in the image by the distance from the pixel in question to the nearest black pixel.

In the synthesized image Is to be processed in this embodiment, a pixel corresponding to the background has a distance value of 0, a pixel adjacent to the background pixel and corresponding to the contour of the object has a distance value of 1, and pixels inside the object have a distance value of 1. As it moves away, the distance value increases. If the pixel value of each pixel in the image is replaced with a distance value and, for example, a pixel with a larger distance value is represented by a higher luminance value, a distance image can be obtained in which the distance from the background is represented by gradation.

FIG. 10 is a diagram showing an example of a distance image, and this distance image Id shows the result of distance transforming the synthesized image Is shown in FIG. Through distance conversion, a distance image Id in which the distance from the background is reflected in the luminance value of each pixel is obtained for the synthetic image Is. The distance image Id is a grayscale image in which the background portion is dark and the luminance increases from the periphery of the object toward the inside. Therefore, an object that has a large area and is nearly circular has a large distance value in the center, and the luminance of the distance image is high (bright). On the other hand, an object with a small area, a thin object, or an object with large unevenness in the surrounding area has a small distance value in the central portion, so the luminance value is low (dark).

The sphere Sp, which is the object to be processed, has a relatively large shape close to a circle in the image. Therefore, there are pixels with relatively large distance values inside the object. On the other hand, for small objects such as floating cells, the distance value for each pixel within the object is smaller. Using this fact, each object in the image is classified into the sphere Sp and other unnecessary objects, and the unnecessary objects are deleted from the range image Id and the synthesized image Is (step S308).

FIG. 11 is a diagram showing the concept of removing unnecessary objects based on distance transform results. As shown in FIG. 11A, the position of each pixel in the distance image Id is represented by coordinates (X, Y), and the distance value of each pixel is plotted on the coordinate plane. In the drawing, the X coordinate and the Y coordinate are represented by one coordinate axis, but more strictly, the distance value axis is perpendicular to the XY coordinate plane.

When expressed in this way, a region with a distance value of 0 corresponds to the background, and one isolated object in the background is represented as a region in which pixels having a distance value greater than 0 are continuous. In the example of FIG. 11A, two objects Obj1 and Obj2 are isolated from each other. The object Obj1, which is relatively large in size and has a uniform shape with few unevenness, has a large peak of the distance value, and the smaller, thin object Obj2 has a small peak. For this reason, an appropriate threshold Vth is set for the distance value, and the object Obj1 whose peak distance value is greater than the threshold Vth is regarded as a sphere and left in the image. Obj2 is regarded as an unnecessary object and deleted from the image.

In this way, relatively large objects, nearly circular objects, etc. remain in the image, while small objects, thin objects, and objects with sharply uneven surroundings are erased from the image. FIG. 11(b) shows an example of the distance image Ida after removing unnecessary objects. At the positions labeled A to D, the small objects and elongated objects that appeared in FIG. 10 do not appear. That is, such objects have been erased. In order to enable accurate unnecessary object removal according to the state of the image, it is preferable that the threshold value Vth is not set to a fixed value, but set to be changeable as a user setting parameter, for example.

In this way, a distance image from which unnecessary objects are excluded and a central image Ica corresponding thereto can be created (step S309). The central image Ica here is an image representing the position of a pixel having a maximum distance value and having a maximum value larger than a predetermined value in the distance image Ida.

FIG. 12 is a diagram showing an example of a central image. More specifically, FIG. 12(a) is the central image Ica corresponding to the distance image Ida shown in FIG. 11(b). 12(b) is a reference diagram in which the central image Ica of FIG. 12(a) is superimposed on the synthesized image Is shown in FIG. In addition, in FIG. 12(b), in order to enhance the visibility of each point indicated by the central image Ica, the position is indicated by a black circle mark.

As can be seen from these figures, and as exemplified with symbols A and B in FIG. be. In this sense, the central image can be said to represent the central position of each object. Therefore, a point appearing in the central image Ica is called a "central point". Although the image in FIG. 12(b) also includes objects that do not include the center point, these are considered unnecessary objects and are excluded from targets for deriving the center point.

As exemplified by symbols C and D in FIG. An object has multiple center points. In other words, an object having multiple center points is considered to be a combination of multiple objects corresponding to the number of center points.

Therefore, an object having a plurality of center points is divided into objects each having one center point by setting appropriate boundary lines between the center points. By doing this, you can split a combined object. Such a concept is nothing but an image processing for object segmentation called "watershed method". That is, based on the information of the distance image Ida and the central image Ica obtained by distance transformation and elimination of unnecessary objects, object division is performed on the composite image Is by the watershed method (step S310), and the images are combined within the image. Object is split.

FIG. 13 is a diagram showing the concept of object division. As schematically shown in FIG. 13(a), an object Obj3 having two center points C1 and C2 is divided into two by setting a boundary line illustrated by a dashed line between the center points C1 and C2. divided into One image processing algorithm for properly setting this boundary is the watershed method. The same can be done for objects containing more than two center points.

FIG. 13(b) is an example of an image after object division. 9 is an image Io showing the outline of an object obtained by dividing the composite image Is shown in FIG. 9 into objects. However, it is assumed that unnecessary objects have been previously deleted from the synthesized image Is. As exemplified by symbols A to C, a plurality of objects that were combined in the synthesized image Is are separated, and as exemplified by symbol D, an object having a conspicuous constriction is also separated from its constricted portion. is divided into two from Each object thus obtained can be regarded as one sphere.

An index is created for the spheres Sp extracted from one frame image If in this way (step S311). For example, it is possible to create an index that associates a unique identification symbol attached to each sphere Sp with information indicating the position of its representative point (for example, center). In addition, the index may include information representing, for example, the area of each sphere Sp, the contour shape, the perimeter, and the like.

In this way, sphere extraction is completed for one frame image If. If there is a frame image to be processed next (YES in step S312), each frame image is sequentially processed by repeating the process from step S301. If the processing is completed up to the last frame image (NO in step S312), the frame extraction processing is completed.

By uniformly executing the processing of steps S301 to S311 for each frame image If, it is possible to extract the sphere Sp in the frame image If. In this case, no association is made between spheres extracted from a plurality of frame images. However, in practice, spheres extracted from two images captured with a small time difference can be associated in a substantially one-to-one relationship. In other words, each of the spheres contained in one image should be present in the other image at approximately the same position and in the same shape, even if there are some changes in position and shape. .

For the purpose of detecting sphere beats, it is necessary to track how the position of a single sphere changes over time. Therefore, it is necessary to know which sphere in one image corresponds to which sphere in the other image. From this point of view, the processing for the images of the second and subsequent frames can be modified as follows.

That is, in the modified sphere extraction process, the process for the image of the second frame omits the generation of the distance image and the center image described as step S309. Then, in object division by the watershed method in step S310, the information of the distance image and the center image generated by the processing of the image of the first frame is used. That is, each object in the image of the second frame is assumed to have its center point at the same position as the object in the image of the first frame. The reason for doing so will be described below.

FIG. 14 is a diagram showing the correspondence of spheres between different frames. FIG. 14(a) shows an example of isolated spheres in the image. It should be noted that FIG. 14(a) shows only a part of the frame image, and the grid of dotted lines is a virtual one added to make it easier to grasp the position of the sphere in the image. The sphere Sp included in the N-th frame image shown on the left side of the figure appears at a slightly displaced position in the (N+1)-th frame image shown on the right side of the figure. Therefore, the center point C(N+1) of the sphere Sp in the (N+1)th frame image may be slightly displaced from the position of the center point C(N) of the sphere Sp in the Nth frame image.

However, the amount of displacement during a short imaging interval is limited, and the center point C(N) in the N-th frame image is contained within, or extremely close to, the sphere Sp in the (N+1)-th frame image. expected to be in close proximity. Therefore, the sphere containing the center point C(N) in the (N+1)th frame image or the sphere closest to the center point C(N) in the (N+1)th frame image is the sphere in the Nth frame image. It can be said that they are identical. This allows correlating corresponding spheres between different frames.

FIG. 14(b) shows an example of spheres combined in an image. Assume that the image of the N-th frame contains two spheres Sp1 and Sp2, which are combined within the image. Let C1(N) and C2(N) be the center points of the spheres Sp1 and Sp2, respectively. Also, the dashed line A represents the sphere boundary defined by the object division at this time.

Assume that the sphere Sp2 is displaced to the right in the image of the (N+1)th frame, as indicated by the dotted line in the drawing. Assuming that the center point of the sphere Sp2 at this time is C2(N+1), it is considered that the center point C2(N+1) is also displaced to the right in response to the movement of the sphere Sp2. On the other hand, the sphere Sp1 hardly moves, and the center point C1(N+1) in the (N+1)th frame is almost the same as the center point C1(N) in the Nth frame.

In the image of the (N+1)th frame, when object division is performed based on the positions of the detected center points C1(N+1) and C2(N+1), as indicated by the dashed line B, the distance between the center points increases. As a result, the boundary is also shifted to the right. As a result, the sphere Sp1 appears to spread rightward. However, the sphere Sp2 actually moved, and the area occupied by the sphere Sp1 remained almost unchanged. In order to suppress the divergence between the object division result and the actual state, the information of the center points C1(N) and C2(N) specified in the Nth frame should be used also in the (N+1)th frame. Just do it.

There are several possible ways of thinking about which frame to use the information of the central image specified in when object segmentation is performed on the images after the (N+2)th frame. For example, it is conceivable to apply the central image information in a single frame, eg the first frame, to the processing of all frames. Also, it is conceivable to apply the information of the center image in the N-th frame imaged immediately before to the processing of the (N+1)-th frame. Furthermore, the central image may be updated every fixed number of frames, and each frame may be processed using the same information until the next update.

In this way, when segmenting an object for one frame image, by executing the watershed method using the distance information and the center information of the frame images captured at different times, the corresponding spheres can be associated between those frame images. It can be carried out. For example, when creating an index for the spheres detected in each frame image, if a common identification symbol is attached to the corresponding spheres, it is possible to clearly record the correspondence of the spheres between images.

As described above, the sphere detection processing of the present embodiment includes the following elemental techniques for detecting spheres from each frame image captured continuously. First, in order to prevent the objects at the edges of the image from being excluded, the peripheral edge of the image is set to a density that can be regarded as the background so that the objects at the edges of the image are surrounded by a closed contour.

In addition, a logical AND synthesized image of a contour region image based on the contour of an object extracted from the original frame image by edge detection and a binarized image obtained by binarizing the frame image based on density (luminance) is an image. indicates the boundary between an object such as a sphere and the background. As a result, the background area surrounded by the object can be correctly recognized as the background.

Then, by eliminating objects that are likely not to be spheres based on the results of the distance transformation and performing object segmentation using the watershed method on the remaining objects, it is possible to identify each of the combined spheres in the image as individual objects. can. At this time, by performing object segmentation using the results of distance transformation in images captured at different times, spheres that move over time can be associated between images, and spheres without motion and spheres with motion can be associated. are combined, they can be divided accurately and stably.

<Thinning process>
FIG. 15 is a diagram showing the principle of thinning processing. As shown in FIG. 15(a), a plurality of original images Ig1, Ig2, . For example, when the frame rate is 60 FPS, the imaging time interval between two original images that are adjacent in time series, that is, that have the shortest imaging time interval is (1/60) second. On the other hand, since the pulsation cycle of the myocardial cell sphere is approximately one second, the displacement of the sphere due to pulsation is extremely small between two adjacent original images. In order to efficiently and reliably detect the displacement of the sphere caused by the pulsation, it is effective to perform image thinning in consideration of the pulsation period.

For example, in order to detect the displacement of a sphere vibrating with a period of 1 second, the imaging time interval between frames is preferably about 1/6 second to 1/3 second. In this case, the original image is thinned out to about 1/10 to 1/20. The degree of this thinning is hereinafter referred to as "thinning rate". In order to flexibly deal with various samples, it is preferable that the thinning rate is not set to a fixed value, but set to be changeable as a user setting parameter, for example. In this case, the original images Ig1, Ig2, .

As shown in FIG. 15(b), thinning is realized by rearranging the images sampled from the original images Ig1, Ig2, . . A group of images If1, If2, If3, . The following processing is executed using each of these images as a substantial frame image. Therefore, in the following explanation of the processing, when simply referring to "frame rate" and "frame image", these refer to those after thinning.

<Beat detection processing>
Next, pulsation detection processing for each sphere detected as described above will be described. The pulsation detection process of this embodiment includes a motion detection step (step S105) and a pulsation determination step (step S106), as shown in FIG. These steps are performed for each sphere in the image.

FIG. 16 is a flow chart showing the motion detection process. Also, FIG. 17 is a diagram showing the relationship between the reference image and the comparison image. In this motion detection step, one image selected from frame images If1, If2, . Step S402), the displacement of the sphere between the two images is obtained, and based on the result, the movement of the sphere between the imaging times of the two images is detected. First, the relationship between the reference image and the comparison image will be described with reference to FIG.

In the motion detection in the frame images If1, If2, . Then, images If2, If3, and If4 picked up after this are sequentially set as comparison images, and the displacement of the sphere between each image and the reference image If1 is obtained.

A comparison between adjacent images in time series reveals the movement of the sphere for a period of time equal to the imaging time interval between frame images. However, the beating period of the spheres is not necessarily constant, and is independent for each sphere. Therefore, the difference in imaging time between the reference image and the comparative image may not be optimal for capturing the motion of the sphere caused by the pulsation. For example, it may happen that the beating sphere happens to be in the same position between the two images. In order to avoid detection omissions caused by this, the comparison images are sequentially replaced, and the displacement of the sphere is detected at several kinds of time intervals. Then, if it is determined that there is motion when the displacement of the sphere is recognized in any combination of the reference image and the comparison image, more reliable motion detection becomes possible.

Here, three comparative images are set for one reference image, but the number of images used as comparative images is not limited to this and is arbitrary. However, if the number of comparison images is increased, the amount of calculation processing also increases and the load on the execution device increases. It is desirable to prevent such omissions in detection.

After applying the three comparison images to the reference image If1, the reference image is then switched to the image If2, and similarly, the three images If3, If4, and If5 are sequentially set as the comparison images to perform the displacement of the sphere. is required. After that, similarly, while switching the reference image to the images If3, If4, . In the flowchart of FIG. 16, this is achieved by loop processing. By doing so, it is possible to detect the movement of the sphere that pulsates at various cycles without omission.

Returning to FIG. 16, the description of the processing contents of the motion detection step is continued. When one reference image and one comparison image are selected (steps S401 and S402), contrast correction is performed to reduce detection errors due to differences in contrast between the two images (step S403). A known correction process can be applied. Also, if correction is not necessary, it can be omitted.

Subsequently, a difference image between the reference image and the comparison image is generated (step S404). The differential image is further binarized (step S405). If there is no motion in the sphere between the two images, there should be no change in image content and the difference should be zero. On the other hand, the larger the displacement, the larger the difference. Therefore, the binarized image of the difference image has information for quantitatively representing the magnitude of displacement.

FIG. 18 is an example of an image showing the difference between the reference image and the comparison image. Consider a reference image I1 and a comparison image I2 containing spheres Sp corresponding to each other. In the figure, the dotted grid is a virtual one added to make it easier to grasp the position of the sphere in the image. Compared to the reference image I1, the sphere Sp in the comparative image I2 is slightly displaced to the right while changing its outer shape. Finding the difference between these yields image I3. In the difference image I3, there appear an approximately annular image caused by the displacement of the outer periphery of the sphere Sp and an image caused by the change in the texture inside the sphere contained therein. An actual difference image is a multi-tone image, and is binarized to facilitate quantitative handling, thereby obtaining a binarized image I4.

A known contraction/expansion process is performed on the binarized image I4 (step S406). As a result, minute speck-like noise components are removed from the binarized image I4. Note that the binarized image after noise removal is hereinafter referred to as a "shrinkage/expansion image" and denoted by reference numeral I4a. This binarized image I4a includes a difference image resulting from the displacement of the outer periphery of the sphere and a difference image resulting from the texture change inside the sphere. Among them, the image caused by the texture change is a factor of error in quantitatively determining the amount of displacement of the sphere. The reason is as follows.

Changes in internal texture are also caused by cellular activity that is not directly related to displacement. Therefore, the sphere may not be displaced and only the internal texture may change. It is necessary to prevent such changes from being mistaken for displacement. Also, even if the internal texture does not change, the entire sphere is displaced, and if the difference between the pixels at the same position is obtained, a significant difference will appear. Since the magnitude of this difference depends on the internal structure of the sphere, the value of the difference is affected not only by the amount of displacement but also by the internal structure of the sphere. More specifically, even with the same amount of displacement, the magnitude of the difference changes depending on the internal structure of the sphere.

Therefore, in order to quantitatively obtain only the displacement of the sphere, it is necessary to eliminate the influence of the internal structure of the sphere and extract only the difference caused by the displacement of the outer periphery of the sphere. For this reason, in this embodiment, a method is adopted in which a search range is set in the vicinity of the outer periphery of the sphere (step S407), and the difference between images within this search range is obtained as a significant value. Two aspects of the search range setting method will be described below, but either method may be applied.

FIG. 19 is a diagram showing a first mode of the search range setting method. More specifically, FIG. 19(a) is a flow chart showing the first aspect of the search range setting method, and FIG. 19(b) is a diagram schematically showing an example of images in the process. As shown in FIG. 19A, contour drawing is performed for each of the sphere in the reference image and the sphere in the comparison image (steps S501 and S502). The contour of the sphere in each frame image has already been obtained in the previous sphere extraction process. Contour drawing is performed based on the contour information obtained in this manner, and the line width is set to a value larger than one pixel. In other words, the outline is drawn in a dilated state.

As a result of the contour drawing, as shown in FIG. 19B, contour images I5 and I6 are obtained by expanding the contours (dotted lines) of the spheres specified in the reference image and the comparative image to a predetermined line width and drawing them. . When these are OR-synthesized (step S503), as shown in the lower part of FIG. A ring-shaped band appears that extends to . This band-shaped area is set as the search range Rs. Here, the image I7 is called a "search range image" in the first mode.

In order to reliably capture the change in the position of the outer peripheral portion as a difference, it is desirable that the line width of the contour to be expanded is thick. On the other hand, if the search range is widened in this way, changes other than changes in the outer periphery of the sphere, such as changes in internal textures in regions close to the outer periphery, will also be detected. Therefore, it is preferable that the line width is not set to a fixed value, but set to be changeable as a user setting parameter, for example.

FIG. 20 is a diagram showing a second aspect of the search range setting method. More specifically, FIG. 20(a) is a flow chart showing the second aspect of the search range setting method, and FIG. 20(b) is a diagram schematically showing an example of images in the process. In order to facilitate understanding of the concept, a simple elliptical shape is used as the outer shape of the sphere in FIG. 20(b). In this method, images I8 and I9 representing the inside of the contour of the sphere are drawn for each of the reference image and the comparison image based on the previously obtained contour information (steps S601 and S602). These images I8 and I9 are subjected to exclusive OR synthesis (step S603).

In the resulting image I10, only the area within the contour of the sphere in only one of the two images is represented in light color (white), and the other areas are represented in dark color (black). When dilation processing is performed to dilate this light-colored region by a predetermined width (step S604), the post-processing image I11 includes both the outlines (dotted lines) of the spheres in the reference image and the comparison image, and the outlines around those outlines. A ring-shaped band appears that extends to . This band-shaped area is set as the search range Rs. Here, the image I11 is called a "search range image" in the second mode.

The search range image representing the search range Rs generated by any one of the search range setting methods is used as a mask to act on the eroded and expanded image I4a. This is equivalent to ANDing the eroded and expanded image I4a and the search range image I7 or search range image I11 (step S408).

FIG. 21 is a diagram showing a logical product synthesized image of an erosion/expansion image and a search range image. Here, the image I7 of the first mode will be used as the search range image. By logically synthesizing the erosion/expansion image I4a representing the difference between the reference image and the comparison image in the outer periphery and inside of the sphere and the search range image I, the synthetic image I12 can be obtained in the vicinity of the outer periphery of the sphere. Only the image corresponding to the change is shown, and the image inside the sphere that existed in the eroded and expanded image I4a is masked and does not appear in the synthesized image. This makes it possible to selectively extract only the position change of the outer periphery, that is, the displacement of the sphere. This composite image I12 is hereinafter referred to as a "motion difference image".

Quantification of displacement is performed as follows. That is, a motion difference area is obtained from the motion difference image I12 (step S409). The motion difference area is obtained as the area of the white portion in the motion difference image I12, and can be simply represented by the number of white pixels in the image. The larger the displacement of the sphere between the reference image and the comparative image, the larger the motion difference area. Therefore, the motion difference area can be used as a numerical value that quantitatively represents the magnitude of the displacement.

However, the motion difference area is affected by the size of the sphere. That is, even if the amount of displacement is the same between a sphere with a large area and a long perimeter and a smaller sphere, the movement difference area that appears as a result of the displacement is naturally larger for the larger sphere. . As an index that can correct this and express the magnitude of displacement regardless of the size of the sphere, the following formula is used:
S = (motion differential area)/(average contour length of sphere)
defines a motion score S. Here, "average edge length of sphere" is the average value of the edge length of the sphere in the reference image and the edge length of the sphere in the comparison image. The contour length of the sphere can be derived from the result of the sphere extraction process, and can be represented, for example, by the number of pixels corresponding to the contour.

The motion score S defined in this way is the difference between images due to the displacement of the sphere normalized by the size of the sphere, and is an index of the magnitude of the displacement regardless of the size of the sphere. .

For each combination of reference and comparison images, a motion score S is calculated for all spheres in the image (step S410). Then, as described above, while sequentially changing the combination of the reference image and the comparison image (steps S411 and S412), the motion score S of the sphere is calculated for each combination. If the motion score S is greater than a predetermined threshold, it can be considered that the sphere exhibits significant motion. Spheres for which movement is detected are preferably recorded along with the direction of movement.

Returning to FIG. 1, next, for each sphere for which the motion score S has been calculated, pulsation determination is performed based on the calculation result (step S106). The presence or absence of pulsation of the sphere can be determined, for example, as follows.

FIG. 22 is a flow chart showing the pulsation determination process. More specifically, FIG. 22 shows processing for determining the presence or absence of pulsation for one sphere. Since each sphere is assigned an index containing a unique identification number, the movement of one sphere can be grasped by extracting the spheres with the same identification number from the detection results in each frame. By arranging the motion detection results for one sphere obtained in the motion detection step in chronological order (step S701), it becomes clear what kind of motion the sphere exhibits over time.

Beating is the continuous and periodic motion of the sphere. Therefore, if the sphere is beating, multiple motions should be detected for the sphere, including those with opposite directions of displacement. Therefore, if the sphere to be determined satisfies these conditions (YES in both steps S702 and S703), it is determined that the sphere is beating (step S704). On the other hand, if the movement is detected only once (NO in step S702) or if the direction of movement is the same a plurality of times (NO in step S703), it is determined that the sphere is not beating. (step S705). In this way, it is possible to determine the presence or absence of pulsation in each sphere. Whether or not the motion directions are opposite can be determined, for example, as follows.

FIG. 23 is a diagram showing an example of determination of the motion direction. When the direction of the detected motion is expressed as a vector, as shown in FIG. 23(a), if the angle θ between the direction of one vector V1 and the direction of the other vector V2 is smaller than 90 degrees, These movements can be judged to be in the same direction. Such movement is not due to pulsation, but simply unidirectional movement of the sphere. On the other hand, as shown in FIG. 23(b), when the angle .theta. The motion direction of the sphere may be obtained from the motion detection result described above, or may be obtained, for example, from the temporal change in the position of the center of gravity of the sphere.

The determination result regarding the pulsation of each sphere is output in an arbitrary manner (step S107). For example, in the sphere contour image Io shown in FIG. 13(b), it is conceivable to perform display output after processing such that the presence or absence of pulsation of each sphere can be identified, for example, by color coding. Further, the frame image If shown in FIG. 2 may be subjected to processing for distinguishing spheres that pulsate, spheres that do not pulsate, and objects other than spheres. Further, detailed information such as presence/absence of pulsation of each sphere, its amplitude, and cycle may be output as a list together with index information. In addition, it is possible to output the determination result in various modes depending on the application.

As described above, in the pulsation detection process of this embodiment, the displacement of the sphere is obtained based on the difference in image content between the reference image and the comparison image captured at different times, and the movement of the sphere is detected from the result. . Displacement due to pulsation appears near the outer periphery of the sphere, while texture change inside the sphere does not necessarily reflect the presence or absence of pulsation. For this reason, in this embodiment, a search range is set in the vicinity of the contour of the extracted sphere, and the image difference within this range is used as an index of the magnitude of displacement. Therefore, it is possible to efficiently obtain the displacement of the sphere due to pulsation without being affected by changes in the internal texture.

In addition, since the pulsation period of the sphere is not necessarily constant and varies among individuals, a plurality of images captured at different times are sequentially applied as comparison images to be combined with the reference image. Therefore, it is possible to reliably detect the pulsation of the sphere that can occur in various cycles.

The magnitude of sphere displacement is evaluated using a motion score S obtained by normalizing the area of the region appearing in the difference image due to the displacement by the contour length of the sphere. Therefore, robust detection and determination are possible with respect to differences in sphere size.

<Device configuration example>
FIG. 24 is a configuration example of a computer device capable of executing the image processing of this embodiment. The computer device 1 has a general configuration as, for example, a personal computer, and includes a CPU (Central Processing Unit) 10, a memory 14, a storage 15, an input device 16, a display unit 17, an interface 18, a disk drive 19, and the like. ing.

The CPU 10 implements the image processing unit 11 as a functional block for executing the above-described image processing in software by executing a control program prepared in advance. Note that dedicated hardware for realizing the image processing unit 11 may be provided. A memory 14 temporarily stores various data generated in the process of calculation by the CPU 10 . The storage 15 stores the control program to be executed by the CPU 10 as well as the image data of the original image, the processed image data, and the like for a long period of time.

The input device 16 is for receiving instruction inputs from an operator, and includes, for example, a mouse and a keyboard. The display unit 17 is, for example, a liquid crystal display having an image display function, and displays various information such as original images, processed images, and messages to the operator. A touch panel in which an input device and a display unit are integrated may be provided.

The interface 18 exchanges various data with external devices via electric communication lines. A disk drive 19 receives an external recording disk 2 on which various data such as image data and control programs are recorded. Image data, control programs, and the like stored in the recording disk 2 are read by the disk drive 19 and stored in the storage 16 . The disk drive 19 may have a function of writing data generated within the computer device 1 onto the recording disk 2 .

The control program for causing the computer device 1 to execute the image processing of the present embodiment may be read out by accessing the recording disk 2 on which it is recorded by the disk drive 19. It may be an aspect provided by a device. The same applies to the original image data.

<Others>
As described above, in the present embodiment, the sphere Sp corresponds to the "cardiomyocyte colony" and the "colony region" of the present invention. Further, the images If1, If2, . The reference image I1 and the comparison image I2 in the pulse detection process are part of the original image and correspond to the "first image" and the "second image" of the present invention, respectively.

The present invention is not limited to the above-described embodiments, and various modifications other than those described above can be made without departing from the spirit of the present invention. For example, in the above embodiments, imaging of the sample, extraction of spheres, and determination of pulsation are described as a series of processes, but such processes are not always required. For example, sphere extraction and pulsation determination may be performed using an image of a sample captured in the past. Also, sphere extraction and beat determination may be performed separately.

Further, in the above embodiment, motion detection is performed by obtaining the displacement of the sphere in the image, but the motion of the sphere due to pulsation may appear as periodic contraction/relaxation. Since the above-described embodiment detects motion from the magnitude of the image difference, it is possible to detect even such contraction and relaxation.

Further, for example, the sphere extraction process in the above embodiment is a process of identifying a relatively circular object such as a sphere from each captured image, and can also be applied to a single still image. is. In addition, the object and purpose of use are not limited to myocardial cell spheres and their pulsation determination, and can be applied to various cell colonies.

Further, for example, when executing the pulsation determination process in the above embodiment, it is necessary to specify the sphere in the image in advance, but the process for specifying the sphere is not limited to the above. The pulsation determination process described above may be performed using sphere information specified by other methods.

Further, for example, in the above-described embodiment, the sphere extraction process is executed for each captured frame image, and the pulsation detection process is executed after thinning out the frames. Alternatively, for example, the captured original image may be thinned out, and the thinned frame image may be used to perform the sphere extraction process and the pulsation detection process.

As described above with reference to specific embodiments, in the image processing method according to the present invention, the binarized image obtained by binarizing the difference image between the first image and the second image, and the search range A value corresponding to the area of the region corresponding to the logical product with the search range image corresponding to may be used as the index value of the magnitude of displacement. With such a configuration, the magnitude of colony displacement can be quantified as a difference, and only the displacement in the vicinity of the contour can be reflected in the index value. A suitable index value can be obtained.

In this case, a value obtained by normalizing the area by the average value of the contour length of the target colony region in the first image and the contour length of the target colony region in the second image may be used as the index value of the displacement. According to such a configuration, even if the amount of displacement is the same, colonies with a larger area have a larger difference and a larger index value, which solves the problem.

For one first image selected from a plurality of original images, a plurality of original images captured at different times are sequentially set as second images one by one, and the amount of displacement of the colony region is detected for each combination. may According to such a configuration, it is possible to avoid omission of pulsation detection even for a sphere whose pulsation period is not necessarily uniform.

Alternatively, a plurality of original images captured at different times may be sequentially set as the first image, and one of the original images captured after the first image may be set as the second image. According to such a configuration, it is possible to detect without omission the movement of the colony region that occurs at any timing during imaging.

Further, the moving direction of the center of gravity of the target colony region between the first image and the second image is defined as the direction of displacement, and displacements larger than a predetermined threshold value are detected multiple times for the target colony region, and those displacements are: It may be determined that a colony corresponding to the target colony region has a pulsation when the directions are opposite to each other. According to such a configuration, it is possible to prevent erroneous determination of movement in one direction, which is different from pulsation, as pulsation. For example, two displacements can be regarded as having opposite directions when the angle formed by the direction vectors indicating the respective displacement directions is greater than 90 degrees.

Further, the step of extracting the colony region from the original image includes a step of extracting a closed curve having a predetermined edge strength from the original image, an image showing the area surrounded by the extracted closed curve, and a binarized image of the original image. and a step of performing object segmentation by the watershed method on objects included in the synthesized image. According to such a configuration, even in a sample in which cell colonies are contained in high density and are in contact with each other, the presence or absence of pulsation can be individually determined for each colony while accurately separating them. is possible.

INDUSTRIAL APPLICABILITY The present invention is applicable to biochemical and medical fields in which cells are imaged and evaluated, and is particularly suitable for evaluation of cardiomyocytes differentiated and cultured from iPS cells or ES cells, for example.

1 computer device 2 recording disk (recording medium)
I1 Original image (first image)
I2 Original image (second image)
I3 difference image I4 binarized image I7, I11 search range image Rs search range Sp sphere (colony, colony area)

Claims (11)

acquiring a plurality of original images of a sample containing colonies of cultured myocardial cells at different imaging times;
a step of extracting a colony region corresponding to the colony from each of a first image and a second image that are different from each other among the original images;
detecting the displacement of the colony region between the first image and the second image, and determining the presence or absence of pulsation of the colony based on the detection result;
for one target colony region among the colony regions, detecting the magnitude of the difference between the first image and the second image within a predetermined search range as the magnitude of the displacement of the target colony region;
The search range includes an annular area obtained by expanding the outline of the target colony area in the first image to a predetermined width, and an annular area obtained by expanding the outline of the target colony area in the second image to the predetermined width. An image processing method, which is an area corresponding to the logical sum of . acquiring a plurality of original images of a sample containing colonies of cultured myocardial cells at different imaging times;
a step of extracting a colony region corresponding to the colony from each of a first image and a second image that are different from each other among the original images;
detecting the displacement of the colony region between the first image and the second image, and determining the presence or absence of pulsation of the colony based on the detection result;
for one target colony region among the colony regions, detecting the magnitude of the difference between the first image and the second image within a predetermined search range as the magnitude of the displacement of the target colony region;
The search range is an annular region including at least a region corresponding to the exclusive OR of the internal region of the target colony region in the first image and the internal region of the target colony region in the second image. Image processing method. A value based on the area of the region corresponding to the logical product of the binarized image obtained by binarizing the difference image between the first image and the second image and the search range image corresponding to the search range is calculated as the displacement. 3. The image processing method according to claim 1, wherein the index value is a magnitude of . A value obtained by normalizing the area with an average value of the contour length of the target colony region in the first image and the contour length of the target colony region in the second image is taken as a value based on the area. 4. The image processing method according to claim 3, wherein the index value of . For one of the first images selected from the plurality of original images, a plurality of the original images having different imaging times are sequentially set as the second image one by one, and the colony region is formed for each combination. 5. The image processing method according to any one of claims 1 to 4, wherein the amount of displacement is detected. 5. A plurality of original images captured at different times are sequentially set as the first image one by one, and one of the original images captured after the first image is set as the second image. The image processing method according to any one of . In the step of determining the presence or absence of pulsation of the colony,
The movement direction of the center of gravity of the target colony region between the first image and the second image is the direction of the displacement,
When displacements larger than a predetermined threshold value are detected multiple times for the target colony region, and the displacements include displacements in opposite directions, the colony corresponding to the target colony region is pulsating. 7. The image processing method according to claim 5 or 6, wherein: 8. The image processing method according to claim 7, wherein two displacements are regarded as opposite directions when an angle formed by directional vectors indicating respective displacement directions is greater than 90 degrees. The step of extracting the colony area from the original image includes:
extracting a closed curve having a predetermined edge strength from the original image;
creating a composite image corresponding to a logical product of an image showing the extracted area surrounded by the closed curve and an image obtained by binarizing the original image;
9. The image processing method according to any one of claims 1 to 8, further comprising the step of performing object segmentation by a watershed method for objects included in said synthesized image. A computer program for causing a computer to execute each step of the image processing method according to any one of claims 1 to 9. A computer-readable recording medium recording the computer program according to claim 10 .

Images