KR20230060888A - Method and apparatus for stitching medical images - Google Patents

Abstract

Disclosed are a medical image stitching method and device. According to an embodiment of the present disclosure, a medical image stitching method includes the steps of: selecting two consecutive images from among several aligned images in an image set to perform stitching; calculating homography between two selected images using an integrated feature detection algorithm based on feature points estimated based on a plurality of different feature detection algorithms; evaluating the accuracy of the calculated homography through a homographic evaluation algorithm; and determining the homography between the two selected images according to evaluation results to perform stitching between the two selected images. Accordingly, accurate stitching which is less dependent on an overlap area can be performed.

Description

Medical image stitching method and apparatus {METHOD AND APPARATUS FOR STITCHING MEDICAL IMAGES}

The present disclosure relates to a medical image stitching method and apparatus for performing registration to generate a panoramic image from a plurality of X-ray images by using a heterogeneous image stitching algorithm.

In general, image stitching is to create one image by matching several consecutive images.

Among image stitching methods, image stitching methods based on feature detection such as SIFT, SURF, and ORB find feature points such as corners and blobs in an image and Feature points may be extracted using surrounding image information. In addition, image stitching may be performed by finding identical (matching) feature points between images and estimating homography (defining a geometric relationship between images) based on this.

In addition, the image-based stitching method may generally perform image stitching by estimating homography using overall image information (eg, similarity of pixel values (intensity) of the entire image).

Although the above-described stitching methods are mainly used when stitching images obtained from a general camera, these methods are often used as they are in stitching medical images such as X-rays.

On the other hand, although the effectiveness of the minimally invasive surgical method has been proven in the treatment of fractures, medical staff rely on C-arm fluoroscopy to perform surgery in order to check bone fragments in real time. Minimally invasive surgery has the advantage of less infection and faster recovery for patients, but it can cause problems for medical staff, such as excessive radiation exposure and an increased possibility of inaccurate fracture reduction.

The C-arm X-ray system is a general imaging device used to observe the condition of a fractured bone in orthopedic surgery. Using the C-arm, bone fragments are aligned in real time during surgery, and the length and angle of the entire bone are measured to check the degree of fracture reduction. At this time, since the field of view of the C-arm is too narrow to visualize the entire bone, the field of view is enlarged by stitching several images using panoramic X-ray images.

That is, X-ray image stitching is a technique of obtaining a panoramic image by combining several X-ray images. In particular, since X-ray images acquired with C-arms in orthopedic surgery have a narrow field-of-view (FOV), a panoramic image with a wide FOV is required to cover the entire region of interest of the patient. In particular, panoramic X-ray images are very important for fracture reduction because they provide important clues needed to evaluate fracture alignment. In addition, in minimally invasive fracture surgery requiring a small incision, panoramic X-ray images are more important because treatment information can be obtained only with 2-dimensional X-ray images.

An X-ray image and an image obtained from a general camera have different video characteristics, such as content, characteristics of equipment, and objects to be photographed.

A precondition for successful inter-image stitching is that overlapping areas between images should exist sufficiently, and a number of matching feature points should exist within the overlapping areas.

However, in the case of X-ray images, a minimum X-ray image is taken due to radiation exposure, and a wide area is taken by minimizing the overlapping area between X-ray images as much as possible. Since shooting is performed, there are few feature point areas that can be used as feature points such as corners and blobs.

Therefore, if the method used when stitching images obtained from a general camera is applied as it is when performing X-ray image stitching, many cases of failure, such as inaccurate stitching, occur.

The foregoing background art is technical information that the inventor has possessed for derivation of the present disclosure or acquired during the derivation process of the present disclosure, and cannot necessarily be referred to as known technology disclosed to the general public prior to filing the present disclosure.

Prior Art 1: Wang L, et al.: Long bone X-ray image stitching using C-arm motion estimation. Informatik aktuell:202, 2009 Prior Art 2: Bai L, Yang JX, Chen XH, Sun YX, Li XY: Medical robotics in bone fracture reduction surgery: a review. Sensors 19, 2019

An object of an embodiment of the present disclosure is to perform matching to generate a plurality of X-ray images into a single panoramic image by utilizing heterogeneous image stitching algorithms, thereby providing accurate and less dependent overlapping regions. stitching is performed.

The purpose of the embodiments of the present disclosure is not limited to the tasks mentioned above, and other objects and advantages of the present disclosure not mentioned above can be understood by the following description and will be more clearly understood by the embodiments of the present disclosure. will be. Furthermore, it will be appreciated that the objects and advantages of the present disclosure may be realized by means of the instrumentalities and combinations indicated in the claims.

A medical image stitching method according to an embodiment of the present disclosure includes the steps of selecting two consecutive images from among a plurality of aligned images of an image set to be stitched, and estimated based on a plurality of different feature detection algorithms. The step of calculating the homography between the two selected images using an integrated feature detection algorithm based on feature points, the step of evaluating the accuracy of the calculated homography through the homographic evaluation algorithm, and the step of evaluating the homography between the two images selected according to the evaluation result. It may include determining the graphics and performing stitching between the two selected images.

An apparatus for stitching medical images according to an embodiment of the present disclosure includes a memory storing at least one command and at least one processor interworking with the memory to execute the command, and when the command is executed by the processor, the processor Select two consecutive images from among several images aligned in the image set to be stitched, and use an integrated feature detection algorithm based on feature points estimated based on a plurality of different feature detection algorithms to determine the relationship between the two selected images. It may be set to calculate homography, evaluate accuracy of the calculated homography through a homographic evaluation algorithm, determine homography between two selected images according to the evaluation result, and perform stitching between the two selected images.

In addition to this, another method for implementing the present disclosure, another system, and a computer readable recording medium storing a computer program for executing the method may be further provided.

Other aspects, features and advantages other than those described above will become apparent from the following drawings, claims and detailed description of the invention.

According to an embodiment of the present disclosure, by performing matching to generate a plurality of X-ray images into a single panoramic image using a heterogeneous image stitching algorithm, accurate stitching that is less dependent on overlapping areas is achieved. can be made to be performed.

In addition, by performing stitching through a series of processes of integrated feature detection, homographic evaluation, and local image-based estimation, stitching accuracy for medical images with small feature points that can be used as feature points such as corners and blobs is improved. can improve

Effects of the present disclosure are not limited to those mentioned above, and other effects not mentioned will be clearly understood by those skilled in the art from the description below.

1 is a diagram schematically illustrating a medical image stitching system according to an embodiment of the present disclosure.
2 is a schematic block diagram illustrating a medical image stitching apparatus according to an embodiment of the present disclosure.
3 is a diagram for explaining an integrated feature detection algorithm according to an embodiment of the present disclosure.
4 is an exemplary diagram for explaining a homography evaluation method according to an embodiment of the present disclosure.
5 is an exemplary diagram for explaining a local image-based estimation method according to an embodiment of the present disclosure.
6 is an exemplary diagram for explaining a method for selecting a blind area area when estimating based on a local image according to an embodiment of the present disclosure.
7 is a flowchart illustrating a method of stitching medical images according to an exemplary embodiment of the present disclosure.
8 is a flowchart illustrating a local image-based estimation method according to an embodiment of the present disclosure.

Advantages and features of the present disclosure, and methods of achieving them, will become clear with reference to the detailed description of embodiments in conjunction with the accompanying drawings.

However, it should be understood that the present disclosure is not limited to the embodiments presented below, but may be implemented in a variety of different forms, and includes all conversions, equivalents, and substitutes included in the spirit and technical scope of the present disclosure. . The embodiments presented below are provided to complete the present disclosure and to fully inform those skilled in the art of the scope of the present disclosure to which the embodiments belong. In describing the present disclosure, if it is determined that a detailed description of related known technologies may obscure the gist of the present disclosure, the detailed description will be omitted.

Terms used in the present disclosure are only used to describe specific embodiments, and are not intended to limit the present disclosure. Singular expressions include plural expressions unless the context clearly dictates otherwise. In the present disclosure, terms such as "comprise" or "having" are intended to indicate that there is a feature, number, step, operation, component, part, or combination thereof described in the specification, but one or more other features It should be understood that the presence or addition of numbers, steps, operations, components, parts, or combinations thereof is not precluded. Terms such as first and second may be used to describe various components, but components should not be limited by the terms. These terms are only used for the purpose of distinguishing one component from another.

Hereinafter, embodiments according to the present disclosure will be described in detail with reference to the accompanying drawings. In the description with reference to the accompanying drawings, the same or corresponding components are assigned the same reference numerals, and overlapping descriptions thereof are omitted. I'm going to do it.

1 is a diagram schematically illustrating a medical image stitching system according to an embodiment of the present disclosure.

Referring to FIG. 1 , a medical image stitching system 1 may include a medical image stitching apparatus 100 , a C-arm 200 , a server 300 and a network 400 .

According to an embodiment, the medical image stitching apparatus 100 is intended to perform stitching of medical images such as X-rays, and heterogeneous ( It is characterized by using a heterogeneous) image stitching algorithm.

In order to perform X-ray image stitching with feature detection, it is essential to extract feature points that are accurately and densely matched within an overlapping region between images. However, since feature points are greatly affected by the properties and size of overlapping regions in successive X-ray images, accuracy and density of matched feature points between images cannot be guaranteed, especially in the case of medical images such as X-rays.

That is, in order to perform successful stitching of the X-ray images obtained from the C-arm 200, an optimal homography for each image must be estimated.

Accordingly, the apparatus 100 for stitching medical images according to an embodiment may perform heterogeneous stitching of X-ray images according to overlapping regions based on homography evaluation.

That is, the medical image stitching apparatus 100 may estimate homography using an integrated feature detection algorithm to sufficiently acquire matching feature points within a limited overlapping area, and evaluate the estimated homography to improve accuracy. You can check. In this case, when it is determined that the estimated homography is not accurate, the medical image stitching apparatus 100 may re-estimate the homography through a local area around the feature point extracted by the integrated feature detection algorithm.

A key part of the X-ray image stitching process of the medical image stitching apparatus 100 is homography to identify deformation between successive images.

In one embodiment, X-ray image stitching may be performed using additional accessories or hardware. For example, a radiolucent ruler for estimating homography may be used. After the radiolucent region is identified in the image, the homography can be calculated by aligning the scales of the radiolucent ruler.

Also, an external plate with a grid of vertical lines, for example, can be used as an additional accessory. Placing an outer plate under the subject causes the grid to be superimposed on the X-ray image. After identifying the intersections of these lines at the boundaries of the image, the X-ray images can be stitched together by matching the image intersections.

Also, for example, an absolute reference panel can be used as a grid of radio-opaque and absolutely positioned markers. Such an absolute reference panel may consist of a cross surrounded by four quadrants with individual tags in binary format, and may allow X-ray image stitching to be performed after each position in the image is estimated by decoding the position markers.

In addition, in one embodiment, instead of using an accessory, X-ray image stitching may be performed using specially designed hardware such as a Camera-Augmented Mobile C arm (CamC) combining a C-arm and a video camera.

Since the homography directly corresponds to the positional change of the X-ray source according to the movement of the C-arm 200, it can be calculated based on video camera data of the estimated movement of the C-arm 200. Therefore, it may be premised that the geometric relationship between the X-ray source of the C-arm 200 and the video camera should be determined in advance through a calibration process. In addition, a visible planar marker pattern may be required to estimate the motion of the C-arm 200 .

On the other hand, another approach of X-ray image stitching is to directly identify deformation using an image matching method that derives homography only from images without using separate accessories or hardware. Image matching methods can be divided into two general categories: direct matching methods and feature point matching methods.

In the case of a direct method, a pixel-by-pixel comparison of two images may be performed for homography. Stitching between successive images can be achieved by quantifying the similarity of images according to sum-of-square differences (SSD), sum of absolute differences, correlation filters, mutual information, and the like.

A common problem of pixel-by-pixel comparison can be caused by brightness variations due to different radiation conditions, and to avoid such non-uniform distribution of pixel values in X-ray images, the histograms of the images can be matched. In addition, image normalization or normalization measures can be used for X-ray image stitching to reduce inaccuracies in similarity calculations due to different brightness levels.

In addition, since image noise also affects the similarity measurement between images, a minimum average correlation energy filter for medical images can be used to remove noise. That is, using the peak to side-lobe ratio (PSLR), which measures the peak sharpness of the correlation plane, X-rays when the image is degraded by noise Image stitching can be done directly.

Meanwhile, the feature detection method may be performed based on feature points extracted from each X-ray image rather than pixel-by-pixel comparison.

For example, after extracting unique feature points of an image using a feature detector such as scale-invariant feature transform (SIFT), matching feature points between images may be identified by performing feature point matching. Feature point matching may be performed by comparing a descriptor (or feature vector) of one feature point with descriptors of all feature points of another image. In an embodiment, an outlier elimination method such as Random Sample Consensus (RANSAC) may be used for more stable feature point matching. Then, the homography can be estimated from the feature points that are stably matched.

In addition, cross-correlation and conditions for detection coherence along with the distance between descriptors for invariably matching feature points extracted by a feature detector such as, for example, Speed-up Robust Feature (SURF) can be suggested An accurate homography may be estimated based on these matching feature points, and a stitched image may be generated through a weighted average fusion process.

In addition, by applying a hybrid method combining SIFT and SURF to improve the performance of feature point matching, for example, it is possible to build a powerful panoramic X-ray image by supplementing the disadvantages of each detection method.

There are several factors to consider when deciding which image matching method is suitable for X-ray image stitching. Although the direct method is simple and efficient, it can estimate only the homography derived from the limited motion of the C-arm (200). Conversely, feature detection methods are time consuming but can achieve more complex stitching without geometric constraints. However, a limitation of the feature detection method is that it is highly dependent on the overlapping area between images. In performing image stitching, the size of the overlapping area for the images should be at least similar and salient features should be included in the overlapping area. Through these prerequisites, homography for X-ray image stitching can be accurately estimated.

However, feature points representing a specific point in an image cannot be accurately defined unless identification areas such as blobs, edges, or corners are sufficiently included in the overlapping area. As a result, image stitching may fail because only a few matching feature points are extracted, and a quantitative area of the overlapping area for successful X-ray stitching cannot be defined. In order to avoid this problem, there may be a method of acquiring many X-ray images to increase the overlapping area.

In one embodiment, a robust X-ray image stitching process may be performed using heterostitching (different types of stitching methods) together with homography evaluation to reduce the dependency on overlapping regions. That is, the medical image stitching apparatus 100 according to an embodiment may perform integrated feature detection to overcome the fact that only a few matched feature points are extracted in the case of a narrow area or unclear overlapping area.

The medical image stitching apparatus 100 according to an exemplary embodiment considers various discriminative properties in an overlapping region of images, for example, SIFT, SURF, BRISK (Binary Robust Invariant Scalable Key Point), and KAZE. can be used

The apparatus 100 for stitching medical images may be used to estimate a homography by combining all feature points by the feature point detection methods or other feature point detection methods in various combinations in an arbitrary number.

Also, since information about the accuracy and density of matching feature points is not derived from the attributes of overlapping regions, the medical image stitching apparatus 100 may perform homography evaluation.

The medical image stitching apparatus 100 may model the movement of the C-arm 200 with respect to a photographing target while obtaining continuous X-ray images to evaluate homography.

The medical image stitching apparatus 100 may replace the homography estimated by the integrated feature detection algorithm with the homography estimated by the local image-based algorithm when it is determined that the homography is inaccurate as a result of the homography evaluation. That is, the medical image stitching apparatus 100 can estimate a more stable homography based on the comparison of a locally limited area around matching feature points located in an overlapping area, so that robust X-ray image stitching can be performed.

Meanwhile, in one embodiment, an embodiment based on the C-arm 200 is described, but it is natural that other types of X-ray imaging devices may also be used.

In the C-arm 200, an X-ray generator (source) is disposed at one end of the C-shaped arm facing each other with a table therebetween, and an image acquisition unit (detector) is disposed at the other end.

In one embodiment, an X-ray image may be obtained using the C-arm 200, and may be referred to as an X-ray image hereinafter.

Meanwhile, in an embodiment, the medical image stitching system 1 may be implemented by the medical image stitching device 100 and/or the server 300, wherein the server 300 may operate the medical image stitching device 100. It may be a server for operation. That is, the server 300 may be a server that performs a process such as calculating matching between images for medical image stitching.

Also, the server 300 may be a database server that provides data for operating the medical image stitching apparatus 100 . In addition, the server 300 may include a web server or an application server. In addition, the server 300 may include a big data server and an AI server required to apply various artificial intelligence algorithms, a calculation server that performs calculations of various algorithms, and the like, and may include the above-described servers or network with these servers. there is.

The network 400 may serve to connect the medical image stitching apparatus 100 and the server 300 and the C-arm 200 in the medical image stitching system 1 . Such a network 400 may be wired networks such as LANs (local area networks), WANs (wide area networks), MANs (metropolitan area networks), ISDNs (integrated service digital networks), wireless LANs, CDMA, Bluetooth, satellite communication However, the scope of the present disclosure is not limited thereto. In addition, the network 400 may transmit and receive information using short-range communication and/or long-distance communication. Here, the short-range communication may include Bluetooth, radio frequency identification (RFID), infrared data association (IrDA), ultra-wideband (UWB), ZigBee, and wireless fidelity (Wi-Fi) technology, and may include long-distance communication. Communications may include code division multiple access (CDMA), frequency division multiple access (FDMA), time division multiple access (TDMA), orthogonal frequency division multiple access (OFDMA), and single carrier frequency division multiple access (SC-FDMA) technologies. can

Network 400 may include connections of network elements such as hubs, bridges, routers, switches, and gateways. Network 400 may include one or more connected networks, such as a multiple network environment, including a public network such as the Internet and a private network such as a secure corporate private network. Access to network 400 may be provided through one or more wired or wireless access networks. Furthermore, the network 400 may support an Internet of Things (IoT) network and/or 5G communication in which information is exchanged and processed between distributed components such as things.

2 is a schematic block diagram illustrating a medical image stitching apparatus according to an embodiment of the present disclosure.

Referring to FIG. 2 , the medical image stitching apparatus 100 may include a communication unit 110 , a user interface 120 , a memory 130 and a processor 140 .

The communication unit 110 may provide a communication interface required to provide a transmission/reception signal between external devices in the form of packet data in conjunction with the network 400 . In addition, the communication unit 110 may be a device including hardware and software necessary for transmitting and receiving signals such as control signals or data signals to and from other network devices through wired or wireless connections. The communication unit 110 may support various things intelligence communication (internet of things (IoT), internet of everything (IoE), internet of small things (IoST), etc.), machine to machine (M2M) communication, vehicle to V2X (vehicle) to everything communication) communication, D2D (device to device) communication, and the like may be supported.

That is, the processor 140 may receive various data or information from an external device connected through the communication unit 110 and may transmit various data or information to the external device. Also, the communication unit 110 may include at least one of a WiFi module, a Bluetooth module, a wireless communication module, and an NFC module.

The user interface 120 controls an algorithm setting process, such as a process of acquiring a medical image to perform stitching of a medical image, setting parameters such as an integrated feature detection algorithm, a homographic evaluation algorithm, and a local image-based homography estimation algorithm. It may include an input interface for inputting user requests and commands to do so.

Also, the user interface 120 may include an output interface for outputting a result performed by the medical image stitching apparatus 100, that is, a stitched panoramic image or the like. That is, the user interface 120 may output results according to user requests and commands. An input interface and an output interface of the user interface 120 may be implemented in the same interface.

The memory 130 may store various information necessary for the operation of the medical image stitching device 100 and/or control (operation) of the server 300 and control software, and may include a volatile or non-volatile recording medium. can do.

The memory 130 is connected to one or more processors 140 and, when executed by the processor 140, causes the processor 140 to control the medical image stitching apparatus 100 and/or the server 300. (cause) You can save codes.

Here, the memory 130 may include magnetic storage media or flash storage media, but the scope of an embodiment is not limited thereto. The memory 130 may include built-in memory and/or external memory, and may include volatile memory such as DRAM, SRAM, or SDRAM, one time programmable ROM (OTPROM), PROM, EPROM, EEPROM, mask ROM, flash ROM, Non-volatile memory such as NAND flash memory, or NOR flash memory, SSD. It may include a compact flash (CF) card, a flash drive such as an SD card, a Micro-SD card, a Mini-SD card, an Xd card, or a memory stick, or a storage device such as an HDD. Also, information related to an algorithm for performing learning according to an embodiment may be stored in the memory 130 . In addition, various information necessary for achieving the purpose of an embodiment may be stored in the memory 130, and the information stored in the memory 130 may be updated as received from a server or an external device or input by a user. may be

The processor 140 may control overall operations of the medical image stitching apparatus 100 . Specifically, the processor 140 is connected to the configuration of the medical image stitching apparatus 100 including the memory 130, and executes at least one command stored in the memory 130 to control the operation of the medical image stitching apparatus 100. can be controlled overall.

Such processor 140 may be implemented in various ways. For example, the processor 140 may include an application specific integrated circuit (ASIC), an embedded processor, a microprocessor, hardware control logic, a hardware finite state machine (FSM), a digital signal processor Processor, DSP) may be implemented as at least one.

The processor 140, as a kind of central processing unit, may control the entire operation of the medical image stitching apparatus 100 by driving control software loaded in the memory 130. The processor 140 may include any type of device capable of processing data. Here, a 'processor' may refer to a data processing device embedded in hardware having a physically structured circuit to perform functions expressed by codes or instructions included in a program, for example.

3 is a diagram for explaining an integrated feature detection algorithm according to an embodiment of the present disclosure, FIG. 4 is an exemplary diagram for explaining a homography evaluation method according to an embodiment of the present disclosure, and FIG. FIG. 6 is an exemplary diagram for explaining a local image-based estimation method according to an embodiment of the present disclosure, and FIG. 6 is an example diagram for explaining a blind area selection method for local image-based estimation according to an embodiment of the present disclosure.

Referring to FIGS. 3 to 6 , the processor 140 of the medical image stitching apparatus 100 will be described.

Integrated feature detection algorithm

In order to increase the number of feature points and the probability of extracting precise matching feature points in a limited overlapping area, it is important to detect a number of distinct feature points in consideration of various attributes of medical images.

Accordingly, referring to FIG. 3 , the processor 140 may combine a plurality of informative feature points calculated by a plurality of different feature detection algorithms for accurate feature point matching.

The processor 140 may estimate feature points calculated by the corresponding feature detection algorithm in the image by applying each feature detection algorithm to two images to be stitched. For example, when the first feature detection algorithm is applied, the first detected feature points of each image may be output.

In addition, the processor 140 may estimate corresponding feature points that are determined to be the same feature points in the two images by comparing descriptors (values formulaically defining feature points) between the estimated feature points. For example, referring to FIG. 3 , based on the first detected feature points of each image, it can be confirmed that first corresponding feature points between images are displayed connected by a yellow line.

Next, the processor 140 may integrate corresponding feature points estimated by a plurality of different feature detection algorithms. In this case, integration may mean generating one feature point set by collecting all estimated corresponding feature points. For example, a set of corresponding feature points may be connected with a yellow line and displayed as an image of integrated corresponding feature points.

Then, the processor 140 estimates the homography using the integrated corresponding feature points.

That is, each feature detection algorithm may extract different unique points for each algorithm using a unique method for each algorithm based on a unique pattern of neighboring pixels.

For example, while SIFT and SURF perform blob-based feature point detection, BRISK and KAZE feature points may be detected more concentrated on boundaries, edges, or corners. In the distribution of detected feature points, BRISK and KAZE appear locally denser than SIFT and SURF.

Therefore, if multiple feature detection is integrated based on each feature detection algorithm, more feature points can be detected and obtained from various local features of an image within a narrow overlapping region compared to single feature detection. That is, in an embodiment, an integrated feature detection algorithm may be applied to detect an integrated feature point, and this integrated feature detection algorithm is implemented as a combination of four feature detection methods, for example, SIFT, SURF, BRISK, and KAZE. It can be. This can be expressed mathematically below. However, in one embodiment, although integrated feature detection is possible by SIFT, SURF, BRISK, and KAZE, it is not limited thereto, and the number of feature point methods and combinations of feature point methods may be variously applied. That is, in an embodiment, integrated feature point detection may be performed by combining an arbitrary number of various feature point extraction methods.

The processor 140 may acquire a plurality of medical images from the C-arm 200. After arranging the acquired medical images in order, it may select at least two consecutive images.

의 i 번째 이미지에서의 통합 특징점은 아래 수학식 1과 같이 나타낼 수 있다.Image set containing selected images

The integrated feature point in the i-th image of can be expressed as in Equation 1 below.

은 4 개의 특징 검출(SIFT, SURF, BRISK, KAZE)에서 얻은 위치 벡터와 기술자 벡터를 의미할 수 있으며, k는 각 특징 검출 방법의 인덱스, j는 각각의 특징 검출 방법의 검출된 특징점의 인덱스를 나타내고,

는 i 번째 이미지에 대한 특징 검출 방법에서 j 번째 검출된 특징점의 위치 벡터와 기술자 벡터 쌍을 나타낸다.here,

may mean position vectors and descriptor vectors obtained from four feature detection (SIFT, SURF, BRISK, KAZE), k is the index of each feature detection method, and j is the index of the detected feature point of each feature detection method. indicate,

represents a position vector and descriptor vector pair of the j-th detected feature point in the feature detection method for the i-th image.

는 다음 수학식 2와 같이 나타낼 수 있다.A feature point located within the FOV of the C-arm (200) is a meaningful feature in feature point matching. However, according to the feature point detection, the boundary of the FOV can be defined in relation to the distinct unique points of the X-ray image. Therefore, in an embodiment, in order to remove feature points beyond the FOV, the FOV radius of the C-arm 200, the feature point location, and the distance between the center of the image may all be used for distance comparison. Refined integrated feature points

can be expressed as in Equation 2 below.

은 각 특징 검출 알고리즘의 정제된 특징점의 인덱스이고,

는 i 번째 이미지의 각 특징 검출 알고리즘의 정제된 특징점의 총 수를 나타내며,

;

는

의 위치 구성 요소,

는 X-ray 이미지의 중심을 나타내고, R은 C-arm(200)의 FOV 반경,

는 j 번째 특징점의 위치와 이미지의 중심 사이의 거리를 나타낸다. 또한,

은 FOV의 경계에서 특징점을 제거하기 위한 반경 임계값 비율을 나타낸다. here

Is the index of the refined feature point of each feature detection algorithm,

Represents the total number of refined feature points of each feature detection algorithm of the ith image,

;

Is

the positional component of

Represents the center of the X-ray image, R is the FOV radius of the C-arm (200),

represents the distance between the location of the j-th feature point and the center of the image. also,

represents the radius threshold ratio for removing feature points at the boundary of the FOV.

를 정의할 수 있다.In an embodiment, the processor 140 may perform feature point matching for each feature detection algorithm based on sum-of-square differences (SSD) after feature point refinement. A matched feature point index set that satisfies Equation 3 below by comparing the descriptor of the m-th feature point of the i-th image with all other feature points of the i+1-th image.

can define

은

의 기술자 구성요소이고, s는 매칭된 특징점들의 인덱스이며,

는 SSD의 임계값을 나타내고,

는

의 e 번째 요소를 나타내며, P는 기술자의 차원(dimension)을 나타내고, m 및 n은 각각 i 번째 및 i+1번째 이미지에서 매칭된 특징점의 인덱스를 나타내며, 수학식 3에서 i 번째와 i+1번째 영상에서 매칭된 특징점의 위치는 수학식 4와 같이, 각각

와

로 표현될 수 있다.here,

silver

is the descriptor component of , s is the index of matched feature points,

represents the threshold of SSD,

Is

represents the e-th element of , P represents the dimension of the descriptor, m and n represent the indices of matched feature points in the i-th and i+1-th images, respectively, and the i-th and i+1 in Equation 3 The positions of the matched feature points in the second image are as shown in Equation 4, respectively.

and

can be expressed as

은

;

의 가로 및 세로 구성 요소를 나타내고,

는 i 번째와 i+1 번째 이미지 사이의 특징 검출 방법 별 매칭된 특징점 인덱스의 총 수를 나타낸다. 그리고

(1)과

(2)는 각각 매칭된 특징점 인덱스의 첫 번째 요소와 두 번째 요소를 나타낸다. here

silver

;

represents the horizontal and vertical components of

represents the total number of matched feature point indices for each feature detection method between the i-th image and the i+1-th image. and

(1) and

(2) represents the first element and the second element of each matched feature point index.

및

로 표시할 수 있으며, 이는 각각

및

의 MSAC 버전을 나타낸다.In one embodiment, in order to remove outliers from matching feature points, for example, an M-estimator sample consensus (MSAC) algorithm may be used. A variant of the RANSAC algorithm, MSAC is an iterative algorithm used to robustly estimate homography from a subset of internal feature points. Using MSAC, better matched feature points can be obtained. The position of the feature point in the continuous image is

and

, which can be expressed as

and

Indicates the MSAC version of .

는 다음 수학식 5와 같이 추정할 수 있다.Based on the location of feature points, the ith homography from the integrated feature detection algorithm

can be estimated as in Equation 5 below.

Here, T represents the transpose of a vector.

Homographic evaluation

In an embodiment, homography is estimated through an integrated feature detection algorithm, but homography uncertainty may remain when the overlap is narrow. Homography may be inaccurate because high density, that is, well-matched feature points cannot be guaranteed. In order to overcome this limitation, in an embodiment, a homographic evaluation algorithm capable of evaluating homography may be applied.

That is, the processor 140 may evaluate the homography estimated through the integrated feature detection algorithm based on the motion of the C-arm 200 with respect to the subject when a series of X-ray images are obtained.

At this time, while obtaining an X-ray image, the subject of the C-arm 200 may be designated as a specific entity. In terms of FOV, the C-arm 200 is very limited compared to a camera. The relative distance between the subject and the detector of the C-arm (200) is limited and the focal length is very long (about a meter). As a result, it can be said that the FOV of the C-arm 200 is essentially very small. Therefore, the motion of the C-arm 200 must follow the shape of the target to cover the entire shape of the object.

Referring to FIG. 4, since the overall shape of the visualized target of the C-arm 200 in general orthopedics is simple and linear, such as the femur, tibia, and spine, the C-arm 200 moves in approximately one direction along the target. motion is expected. These shifts are characterized by patterns that can be described by mathematical formulas when acquiring successive X-ray images for reconstruction of panoramic images. Accordingly, in an embodiment, homography may be evaluated based on knowledge of the motion of the C-arm 200 .

과 이미지 회전 정도의 기준

을 결정할 수 있다. C-arm(200)이 대상체(bones)를 찍을 때 움직임, 즉 C-arm(200)의 FOV에 따른 이동 정도 및 회전 정도를 추정하고 이를

,

에 적용하여 결정할 수 있다. 도 4에 도시된 바와 같이, 예를 들어 대퇴부 뼈라고 가정하면 C-arm(200)은 대퇴부 뼈 전체를 촬영해야 하기 때문에 대퇴부 뼈의 전체의 근사적인 형태(두꺼운 직선)로 이동하게 된다. 따라서 C-arm(200)에서는 회전이 거의 없이 한 방향으로 이동하면서 촬영하게 된다. 이를

,

결정에 반영하여 호모그래피를 평가할 수 있다.That is, in one embodiment, the movement of the C-arm 200 is modeled to determine the degree of horizontal movement used in the self-evaluation index.

and standard for image rotation degree

can determine When the C-arm (200) takes pictures of the target (bones), the motion, that is, the degree of movement and rotation according to the FOV of the C-arm (200) is estimated, and

,

can be determined by applying As shown in FIG. 4, assuming, for example, the femoral bone, the C-arm 200 moves in an approximate shape (thick straight line) of the entire femoral bone because the entire femoral bone must be photographed. Therefore, the C-arm 200 takes pictures while moving in one direction with almost no rotation. this

,

The homography can be evaluated by incorporating it into the decision.

Meanwhile, in one embodiment, the processor 140 may perform the homographic evaluation in the above two steps, such as a) self-evaluation and b) relative evaluation.

와 회전 값

은 수학식 6과 같이 표현될 수 있다. 이러한 기준에 기반하여, i 번째 호모그래피의 자체 평가 인덱스

가 정의될 수 있다.In the self-evaluation, translation and rotation of homography may be considered. In self-evaluation, the horizontal shift value of the ith homography

and rotation value

Can be expressed as in Equation 6. Based on these criteria, the self-assessed index of the ith homography

can be defined.

는 가중치로, 이미지의 너비일 수 있으며,

및

은 각각 너비와 각도의 여백(margin), 즉 이미지의 수평 이동 정도의 기준 및 이미지 회전 정도의 기준을 나타낼 수 있다. here,

is a weight, which may be the width of the image,

and

may indicate a margin of width and angle, that is, a criterion for a degree of horizontal movement of an image and a criterion for a degree of rotation of an image.

Relative evaluation can be evaluated by comparing the current homography with the previous homography. For example, if there are three input images, the homographies estimated from the first and second images can be compared with the estimated homographies from the second and third images.

와 호모그래피 간의 회전 각도 차이

는 결정 요인이다. i 번째 호모그래피에서 상대 평가 인덱스

는, 현재 이미지에서 아래 수학식 7에 의해 결정될 수 있다.In relative evaluation, the horizontal shift difference between homographies

Rotation angle difference between and homography

is the determining factor. Relative evaluation index in the ith homography

, may be determined by Equation 7 below in the current image.

에 의해 호모그래피의 정확성에 대한 최종 평가가 결정될 수 있으며,

은 수학식 8과 같이 나타낼 수 있다.Since the first homography of the image set had no prior homography, only self-evaluation can be performed. From the second homography, relative evaluation may be performed together with self-evaluation. Through these two evaluation steps, the evaluation index of the ith homography

The final evaluation of the accuracy of the homography can be determined by

can be expressed as in Equation 8.

는 논리적 OR 연산을 나타낸다. 평가 인덱스가 1일 때 호모그래피가 정확한 것으로 판단할 수 있으며, 그렇지 않으면, 호모그래피를 재 추정하기 위해 로컬, 즉 대응하는 이미지 영역을 기반으로 하는 로컬 영역 추정 알고리즘이 사용될 수 있다. 다시 말하면, 프로세서(140)는 자체 평가 인덱스와 상대 평가 인덱스 중에 1이 있으면 호모그래피가 정확한 것으로 판단할 수 있고, 그렇지 않으면 호모그래피가 정확하지 않다고 판단할 수 있다.here,

represents a logical OR operation. When the evaluation index is 1, it can be determined that the homography is correct. Otherwise, a local area estimation algorithm based on a local, corresponding image area, can be used to re-estimate the homography. In other words, the processor 140 may determine that the homography is accurate if 1 is present between the self-evaluation index and the relative evaluation index, and otherwise may determine that the homography is not accurate.

Local image-based estimation

In one embodiment, if it is determined that the homography derived through the integrated feature detection algorithm is inaccurate, the homography may be re-estimated by replacing pixel values around the matched feature point (point).

In one embodiment, based on the matched feature points in the overlapping area from the integrated feature detection algorithm, a corresponding local area may be determined as a square area centered on the location of the feature point on each image. Then, the processor 140 may re-estimate the homography based on the peak value of the flatness-based extended phase correlation of the corresponding local area.

At this time, in re-estimating the homography, precisely matched feature points are not required, so the dependence on the overlapping region is greatly reduced, unlike feature point detection.

In one embodiment, local image-based estimation may consist of two parts: corresponding local region selection and flatness-based extended phase correlation.

Local Image Based Estimation - Corresponding local region selection

In one embodiment, matched feature points may be modified in advance to determine the corresponding local area. First, among the matching feature points of the integrated feature detection result, a matched feature point candidate for a local area may be selected in consideration of the width of each image.

Referring to FIG. 5(b), corresponding feature points derived through an integrated feature detection algorithm are shown. Matched points mean corresponding feature points that exist within 3/4 of the width of the first image and 1/4 of the width of the second image among all corresponding feature points derived through the integrated feature detection algorithm. At this time, it appears similar to the red area of FIG. 5(b), which may mean an overlapping area. Also, since the overlapping area is narrow and the exact width of the overlapping area cannot be determined, it can be assumed that the location of the matched feature map is limited to within 1/4 of the image width.

In one embodiment, in selecting a corresponding local area, a rectangular area is set based on the matched point, and the setting of the rectangular area is determined by the position of a feature point (a red circle and a green cross in FIG. 5(b)). marked) as the center of the rectangle, and can be experimentally set to have a length of 5% of the total image size. For example, if the width of the image is 600 pixels, the user can define the size of the rectangle as 30 pixels.

및

는, 수학식 9와 같이 수학적으로 나타낼 수 있다.Location of Matched Feature Points

and

Can be expressed mathematically as shown in Equation 9.

는

의 수평 위치 성분을 나타낸다. 그리고 로컬 영역은 매칭된 특징점 후보를 중심으로 하는 정사각형 영역으로 정의되어, 매칭된 특징점이 중첩 영역에 위치할 수 있다. here

Is

represents the horizontal position component of Also, since the local area is defined as a square area centered on the matched feature point candidate, the matched feature point may be located in an overlapping area.

Referring to FIG. 6, the selection process of the local area is described. The feature maps matched in the integrated feature detection result are included in the two side-by-side X-ray images. The following image shows one matched feature point candidate that satisfies Equation 9, and the corresponding local area is indicated by a red dotted line rectangle in each image. After all local areas satisfying Equation 9 are determined, the determined local areas may become candidates used to calculate the flatness-based extended phase correlation.

Extended phase correlation with flatness based on local image-based estimated flatness

In one embodiment, the homography of local image-based estimation may be approximated to a rigid transformation by reflecting the motion of the C-arm 200 . Therefore, components for translation and rotation are needed.

In one embodiment, after calculating the rotation in the frequency domain, it is possible to calculate the extended phase correlation enabling the translation calculation with sub-pixel accuracy. Rotation can be estimated using the local domain Fourier transform amplitude.

For example, when there are consecutive images, the processor 140 may perform matching rectangles existing within 3/4 of the width of the first image and 1/4 of the width of the second image among all matching rectangle local images. You can only select local images of .

Also, the processor 140 may estimate a rotation component between the first image and the second image using Fourier transform.

Next, the processor 140 may use the estimated rotation value to correct a rotation component that is different between the second image and the third image among the three consecutive images. Accordingly, only a difference in moving components exists between the two images. Accordingly, the processor 140 may estimate motion components of the second image and the third image through phase correlation with flatness.

That is, the processor 140 may perform homography transformation based on the estimated motion component and rotation component.

At this time, the processor 140 performs rotation estimation using the Fourier transform of the image in order to estimate the rotation component. Since the Fourier transform is invariant to the motion component of the image, it is possible to estimate a rotation value by slightly rotating the Fourier transformed image so that they match, and this can be defined as a homography rotation component value of the two images.

In addition, the processor 140 uses flatness-based phase correlation to estimate the motion component. First, the rotation between the two images is matched using the rotation component using the Fourier transform. As a result, only motion components exist between the two images, and motion can be estimated by applying flatness-based phase correlation.

및

를 중심으로 하는 i 번째 및 i+1 번째 이미지의 해당 로컬 영역

및

를 정의하고, 이들의 푸리에 변환은 아래 수학식 10과 같이 나타낼 수 있다.In one embodiment, each

and

Corresponding local areas of the i-th and i+1-th images centered at

and

, and their Fourier transform can be expressed as in Equation 10 below.

는 푸리에 변환을 나타낸다. 이동은 푸리에 변환의 위상에만 영향을 미치므로 이미지의 이동 및 회전 구성 요소를 분리할 수 있다. 따라서 최대 상관 관계와 푸리에 변환의 진폭을 일치시켜 회전을 추정할 수 있다. here

represents the Fourier transform. Since translation only affects the phase of the Fourier transform, we can separate the translation and rotation components of the image. Therefore, rotation can be estimated by matching the maximum correlation with the amplitude of the Fourier transform.

는

와 같은 부등식을 따른다고 가정할 수 있으며, 각도 에서 직사각형 그리드 값의 평균으로 계산될 수 있다. 다음 수학식 11을 참조하면, 추정된 회전 각도

는 상관 관계가 최대화되는 값으로 계산될 수 있다.Also, in one embodiment, the correlation can be converted into polar coordinates to efficiently calculate the correlation, and the amplitude as a function of angle based on the integration along the radial line can be calculated. thus,

Is

It can be assumed to follow the same inequality as , and can be calculated as the average of the rectangular grid values at angle . Referring to Equation 11 below, the estimated rotation angle

can be calculated as a value that maximizes the correlation.

는 극좌표 각도를 나타내며, 또한

는 상관 관계를 나타낼 수 있다.where r is the radius,

denotes the polar coordinate angle, and also

may indicate a correlation.

In an embodiment, flatness-based extended phase correlation may be applied to estimate translation. It is based on the property that the inverse Fourier transform of the normalized cross-power spectrum is the impulse function at the point of the translated pixel value.

및

의 피크 값

은 다음 수학식 12와 같이 계산될 수 있다.In the case of extended phase correlation for estimating subpixel translation, the subpixel translation image can be assumed to be a downsampled image converted by integer pixels, and a closed-form solution of extended phase correlation. ) can be approximated by a two-dimensional sinc function to obtain In addition, the shift to sub-pixel accuracy between corresponding local areas can be indicated by the maximum peak (z-axis) of the sinc function. that local area

and

peak value of

Can be calculated as in Equation 12 below.

및

는 각각 수평 및 수직 방향의 다운샘플링 인수이다.

및

는 이미지를 다운샘플링하기 전의 수평 및 수직 방향의 정수 이동 값이다. 일 실시 예에서는, 이동을 추정할 때 로컬 영역의 평탄도도 고려한다. 플랫(flat)한 로컬 영역의 확장된 위상 상관은 sinc 함수의 최대 피크에 대한 모호성을 유발할 수 있다. 이에, 다음 수학식 13을 참조하여, 이동 추정에서 상기 특징점에 의해 결정되는 로컬 영역을 제외하기 위해, 해당 로컬 영역의 픽셀 값(intensity)

및

의 최소 표준 편차를 정의하는 로컬 영역의 평탄도

를 계산할 수 있다.here

and

are the downsampling factors in the horizontal and vertical directions, respectively.

and

is an integer shift value in the horizontal and vertical directions before downsampling the image. In one embodiment, the flatness of the local area is also taken into account when estimating movement. The extended phase correlation of the flat local region can cause ambiguity about the maximum peak of the sinc function. Therefore, referring to Equation 13 below, in order to exclude the local area determined by the feature point from the motion estimation, the pixel value (intensity) of the corresponding local area

and

The flatness of the local area defining the minimum standard deviation of

can be calculated.

Here, i means the index of the image, min() indicates the minimum value between the two values, and std() indicates the standard deviation. That is, flatness can be defined as a standard deviation value calculated by calculating the standard deviation of pixel values of pixels within the local region of the matching rectangle.

를 계산할 수 있다.In an embodiment, an optimal local region for estimating homography may be determined after calculating peak values and flatness of all corresponding regions in matching feature point candidates. Referring to Equation 14 below, the criterion for the peak value and flatness may be defined, and referring to Equation 15, parallel conversion between consecutive images may be performed.

can be calculated.

는 해당 로컬 영역 간의 추정된 이동을 나타낸다.where p is the location of the matched feature point in the first image between two consecutive images,

denotes the estimated movement between corresponding local areas.

Referring to FIG. 6 , the rightmost rectangular local area means the final rectangular local area. That is, rotation estimation and correction are performed for a plurality of selected matching local rectangular regions, flatness-based phase correlation is calculated, and one local that satisfies Equation 14 (maximum peak value is the largest and satisfies the flatness condition) A region is selected, and through this, the calculated rotation value and movement value can be finally set as homography.

는 수학식 16과 같이 계산될 수 있다. i-th homography in local image-based estimation, using rotations and translations estimated in the local corresponding image-based method

Can be calculated as in Equation 16.

Image blending

In the case of a panoramic image, each X-ray image may be mixed using the estimated homography. After registering images using homography, overlapping regions between consecutive images can be blended into one image. In one embodiment, a simple and fast method may be used, the maximum blending value. After selecting the maximum pixel values from the two images of the overlapping area, the final stitched X-ray image can be obtained.

In summary, in one embodiment, in the selection of a feature detection algorithm for integrated feature detection, a combination of four feature detections, SIFT, SURF, BRISK, and KAZE, is performed to obtain a sufficient number of feature points within the overlapping region of the image. It was explained with an example. As described above, it is not limited thereto, and the number of feature point methods and combinations of feature point methods may be variously applied.

Compared with the image properties of ordinary digital cameras, X-ray images are simple and textureless because the background is empty, so that few distinct features are observed in the background of the overlapping area. Since the C-arm 200 is a medical device used for detailed diagnosis of a subject, an X-ray image is acquired only for a local part of the subject. Therefore, in the X-ray image, a limited number of feature points are detected from single feature detection.

In order to overcome the above problems, in an embodiment, additional feature points reflecting various distinct regions of interest may be detected.

Accordingly, in one embodiment, various unique feature points such as blobs, edges, and boundaries may be detected by applying an integrated feature detection algorithm.

In addition, in an embodiment, if it is determined that the homography estimated based on the integrated feature detection is incorrect after the homographic evaluation, the homography may be re-estimated based on the local region. In an embodiment, the location of the local region may be designated by a matched feature point through an integrated feature detection algorithm in order to re-estimate the homography based on the local correspondence region. After selecting a matched feature point in the overlapping area from among all matched feature point candidates, the local correspondence area may be determined as a square area centered on the location of the selected feature point, and its size may not be fixed.

In an embodiment, the width of the local area may be limited according to the position of the selected feature point. When the position of the selected feature point is near the boundary of the image, the maximum width of the local region can be determined by the distance between the position of the feature point and the boundary of the image. The local area's minimum width can be more flexible than its maximum width.

However, if the size of the local region is too small, the ambiguity of the peak value of the extended phase correlation increases due to the small number of adjacent pixels, and thus the homography may be erroneously estimated. In an embodiment, the width of the local region based on the above factors can be set to a range of 5% or more of the width of the image.

Hereinafter, a series of processes for performing medical image stitching by the processor 140 will be described with reference to FIGS. 7 and 8 .

7 is a flowchart illustrating a method of stitching medical images according to an exemplary embodiment of the present disclosure. Referring to FIG. 7 , in step S100, the processor 140 selects two consecutive images from among several aligned images of an image set to be stitched.

For example, after arranging the input images in order, an image and its successive X-ray images can be selected. In addition, after arranging the images in order for the selected image, the image and its successive X-ray images can be selected.

In step S200, the processor 140 calculates a homography between the two selected images using an integrated feature detection algorithm.

At this time, the processor 140 applies a first feature detection algorithm and a second feature detection algorithm to the selected two images, compares descriptors between feature points of each of the two images estimated through the first feature detection algorithm, and obtains a corresponding first The feature points may be estimated, and corresponding second feature points may be estimated by comparing descriptors between feature points of each of the two images estimated through the second feature detection algorithm.

Also, the processor 140 may integrate the first feature points and the second feature points to extract an integrated feature point set, and calculate a homography using the integrated feature point set.

Also, the processor 140 may remove feature points outside the FOV range of the image capture device based on the FOV radius of the image capture device, the location of the feature point, and the center location of each image. The image capturing device may include the C-arm 200.

That is, in an embodiment, homography may be estimated using integrated feature detection for the selected image. After the distinctive feature points of the image are extracted, homography can be estimated based on the matched feature points, and the estimated homography can be verified through a homographic evaluation method as will be described later.

In step S300, the processor 140 evaluates the accuracy of the homography calculated based on the integrated feature detection algorithm through a homographic evaluation algorithm.

In this case, the processor 140 may perform self-evaluation based on translation and rotation components of the estimated homography between the two images.

To perform self-evaluation, the processor 140 may compare the horizontal shift value in the estimated homography with the weighted image horizontal shift reference value, and compare the rotation value in the estimated homography with the image rotation reference value. there is.

In this case, the processor 140 determines that the homography is correct if the horizontal movement value in the estimated homography is less than the weighted image horizontal movement reference value or the rotation value in the estimated homography exceeds the image rotation reference value. It is determined as not (0), and if not, it may be determined that the homography is correct (1). In an embodiment, the image horizontal movement reference value and the image rotation reference value may be set by modeling the motion of the image capture device.

Also, the processor 140 may perform relative evaluation by comparing the estimated homography between the two images with the estimated previous homography between the previous image and the first image of the two images. That is, in an embodiment, only self-evaluation may be performed when there is no previous image.

The processor 140 may perform relative evaluation based on a difference in horizontal translation values and a difference in rotation values between the estimated homography between the two images and the estimated previous homography between the previous image and the first image of the two images. .

At this time, the processor 140 subtracts the estimated horizontal shift value of the previous homography between the previous image and the first of the two images from the estimated horizontal shift value of the homography between the two images and exceeds 0, and the two images If the value obtained by subtracting the rotation value of the estimated previous homography between the previous image and the first of the two images from the rotation value of the estimated homography between the two images is less than the image rotation criterion value, the homography is correct(1) , and if not, it may be determined that the homography is not correct (0).

Further, the processor 140 may perform final homography evaluation based on the results of self-evaluation and relative evaluation.

The processor 140 determines that the homography estimated by the integrated feature detection algorithm is correct when at least one of the self-evaluation result and the relative evaluation result is output as correct (1), and otherwise In this case, it may be determined that the homography estimated by the integrated feature detection algorithm is not accurate.

In step S400, the processor 140 determines that the evaluation result through the homographic evaluation algorithm and the homography calculated based on the integrated feature detection algorithm are correct (yes in step S300), the calculation based on the integrated feature detection algorithm One homography is determined as the final homography.

That is, in an embodiment, if it is determined that the homography is correct as a result of verifying the homography based on integrated feature detection, the homography based on integrated feature detection may be regarded as an optimal homography for the corresponding image and stored.

On the other hand, as will be described later, integrated feature detection is replaced with local, corresponding image-based estimation, and a homography estimated by applying a local area estimation algorithm is considered an optimal homography and can be stored. This process can be run repeatedly until the last image in the image set has been processed.

That is, in step S500, the processor 140 uses a local region estimation algorithm when it is determined that the homography calculated based on the homographic evaluation algorithm as a result of the evaluation through the homographic evaluation algorithm is not correct (No in step S300). to re-estimate the homography between the two images.

At this time, the processor 140 may select a corresponding local area between the two images in each of the two images, and re-estimate the homography based on the selected local area.

The processor 140 may select, as a local area, a rectangle having a length corresponding to a predetermined ratio of the size of the corresponding image, with the location of the corresponding feature point between the two images as a center in each of the two images.

Also, the processor 140 may re-select and select a local area included in certain areas of the two images among the selected local areas. Based on the selected local region, the rotation component between the two images is estimated, the motion component between the two images is estimated by applying the estimated rotation component, and the homography can be re-estimated based on the estimated rotation and motion components.

In step S600, the processor 140 determines the re-estimated homography based on the local area estimation algorithm as the final homography, and in step S700, stitching between the two images is performed through the determined final homography.

At this time, in the case of a panoramic image, since each X-ray image can be mixed using the estimated homography, the processor 140 registers the image using the finally determined homography, and then divides the overlapping area between consecutive images into one. Images can be blended.

In one embodiment, a simple and fast method may be used, the maximum blending value. After selecting the maximum pixel values from the two images of the overlapping area, the final stitched X-ray image can be obtained.

8 is a flowchart illustrating a local image-based estimation method according to an embodiment of the present disclosure. Referring to FIG. 8, a method of re-estimating homography based on a local area estimation algorithm will be described in more detail.

Referring to FIG. 8 , in step S510, the processor 140 selects a corresponding local area between the two images in each of the two images.

In an embodiment, if it is determined that the homography derived through the integrated feature detection algorithm is inaccurate, the processor 140 may re-estimate the homography by substituting pixel values around matched feature points. Therefore, the processor 140 may determine a corresponding local area as a square area centered on the location of the feature point on each image, based on the matched feature point in the overlapping area from the integrated feature detection algorithm.

In step S520, the processor 140 selects a local area included in a predetermined area of the two images from among the selected local areas.

The processor 140 may modify the matched feature points in advance to determine the corresponding local area. First, among the matching feature points of the integrated feature detection result, a matched feature point candidate for the local area may be selected in consideration of the width of each image. there is.

The processor 140 may select, as a local area, a rectangle having a length corresponding to a predetermined ratio of the size of the corresponding image, with the location of the corresponding feature point between the two images as a center in each of the two images. That is, the processor 140 may reselect and select the local area included in the predetermined area of the two images among the selected local areas.

In step S530, the processor 140 estimates a rotation component between the two images based on the selected local region.

At this time, the processor 140 may estimate a rotation value by applying a Fourier transform to the selected local region and rotating the local region to which the Fourier transform is applied so that the local regions of the two images match.

In one embodiment, the homography of local image-based estimation may be approximated to a rigid transformation by reflecting the motion of the C-arm 200 . Therefore, components for translation and rotation are needed.

Accordingly, the processor 140 may calculate the rotation in the frequency domain and then calculate the extended phase correlation enabling the movement calculation with sub-pixel accuracy. Rotation can be estimated using the local domain Fourier transform amplitude.

In step S540, the processor 140 corrects and matches the second image of the two images through the estimated rotation component, and in step S550, moves between the two images based on phase correlation with flatness. estimate the ingredients.

At this time, the processor 140 may perform rotational component estimation and correction reflecting rotational components for all corresponding local regions of the two images. In addition, the processor 140 applies flatness-based phase correlation analysis to all corresponding local areas of the two corrected images, and selects one corresponding local area that has the largest maximum peak value and satisfies a preset flatness condition. can be derived Here, flatness may be calculated as a smaller standard deviation value among standard deviations of pixel values (intensities) of all corresponding pixels in the local area.

That is, the processor 140 performs rotation estimation using the Fourier transform of the image in order to estimate the rotation component. Since the Fourier transform is invariant to the motion component of the image, it is possible to estimate a rotation value by slightly rotating the Fourier transformed image so that they match, and this can be defined as a homography rotation component value of the two images.

In addition, the processor 140 uses flatness-based phase correlation to estimate the motion component. First, as a result of matching the rotation between the two images using the rotation component using the Fourier transform, only the motion component exists between the two images. and the motion component can be estimated by applying the flatness-based phase correlation method.

In step S560, the processor 140 re-estimates the homography based on the estimated rotation and motion components.

That is, an embodiment is based on heterogeneous stitching that depends on overlapping regions in X-ray images, and when integrated feature detection, homographic evaluation, and local image-based estimation are used, X-ray image stitching performance is less in overlapping regions. Depending on the dependencies, more powerful panoramic C-arm X-ray images can be constructed.

Embodiments according to the present disclosure described above may be implemented in the form of a computer program that can be executed on a computer through various components, and such a computer program may be recorded on a computer-readable medium. At this time, the medium is a magnetic medium such as a hard disk, a floppy disk and a magnetic tape, an optical recording medium such as a CD-ROM and a DVD, a magneto-optical medium such as a floptical disk, and a ROM hardware devices specially configured to store and execute program instructions, such as RAM, flash memory, and the like.

Meanwhile, the computer program may be specially designed and configured for the purpose of the present disclosure, or may be known and available to those skilled in the art in the field of computer software. An example of a computer program may include not only machine language code generated by a compiler but also high-level language code that can be executed by a computer using an interpreter or the like.

In the specification of the present disclosure (particularly in the claims), the use of the term "above" and similar indicating terms may correspond to both singular and plural. In addition, when a range is described in the present disclosure, as including the invention to which individual values belonging to the range are applied (unless otherwise stated), each individual value constituting the range is described in the detailed description of the invention Same as

Unless an order is explicitly stated or stated to the contrary for steps comprising a method according to the present disclosure, the steps may be performed in any suitable order. The present disclosure is not necessarily limited to the order of description of the steps. The use of all examples or exemplary terms (eg, etc.) in this disclosure is simply to explain the present disclosure in detail, and the scope of the present disclosure is limited due to the examples or exemplary terms unless limited by the claims. it is not going to be In addition, those skilled in the art can appreciate that various modifications, combinations and changes can be made according to design conditions and factors within the scope of the appended claims or equivalents thereof.

Therefore, the spirit of the present disclosure should not be limited to the above-described embodiments, and not only the claims to be described later, but also all ranges equivalent to or equivalent to these claims are within the scope of the spirit of the present disclosure. will be said to belong to

1: Medical image stitching system
100: medical image stitching device
110: Communication Department
120: user interface
130: memory
140: processor
200: C-arm X-ray imaging device
300: server
400: network

Claims (20)

A method of stitching medical images in which at least part of each step is performed by a processor,
selecting two contiguous images from among a plurality of aligned images of an image set to be stitched;
calculating homography between the two selected images using an integrated feature detection algorithm based on feature points estimated based on a plurality of different feature detection algorithms;
Evaluating the calculated accuracy of the homography based on a homographic evaluation algorithm; and
Determining homography between the two selected images according to an evaluation result based on the homographic evaluation algorithm and performing stitching between the two selected images,
Methods for stitching medical images.
According to claim 1,
The step of calculating the homography,
applying a first feature detection algorithm and a second feature detection algorithm to the two selected images;
estimating corresponding first feature points by comparing descriptors between feature points of each of the two images estimated through the first feature detection algorithm;
estimating corresponding second feature points by comparing descriptors between feature points of each of the two images estimated through the second feature detection algorithm;
extracting an integrated feature point set by integrating the first feature points and the second feature points; and
Comprising the step of calculating a homography using the integrated feature point set,
Methods for stitching medical images.
According to claim 2,
The step of calculating the homography,
Further comprising removing feature points outside the FOV range of the diagnostic imaging device based on a field of view (FOV) radius of the diagnostic imaging device, a location of the feature point, and a center location of each image,
Methods for stitching medical images.
According to claim 1,
The evaluation step is
performing self-evaluation based on translation and rotation components of the estimated homography between the two images;
performing relative evaluation by comparing the estimated homography between the two images with the estimated previous homography between a previous image and a first image of the two images; and
Performing a final homography evaluation based on the self-evaluation and the relative evaluation results,
Methods for stitching medical images.
According to claim 4,
The step of performing the self-assessment,
comparing a horizontal shift value in the estimated homography with a weighted image horizontal shift reference value; and
Comparing a rotation value in the estimated homography with an image rotation reference value,
Methods for stitching medical images.
According to claim 5,
The step of performing the self-assessment,
If the horizontal shift value in the estimated homography is less than the weighted image horizontal shift reference value, or if the rotation value in the estimated homography exceeds the image rotation reference value, the homography is not accurate (0) , and if the horizontal movement value in the estimated homography is greater than or equal to the weighted image horizontal movement reference value and the rotation value in the estimated homography is less than the image rotation reference value, the homography is correct (1) Including the step of judging,
Methods for stitching medical images.
According to claim 5,
The image horizontal movement reference value and the image rotation reference value are set by modeling the motion of the image capturing device.
Methods for stitching medical images.
According to claim 4,
The step of performing the relative evaluation,
Performing the relative evaluation based on a difference in horizontal movement value and a difference in rotation value between the estimated homography between the two images and the estimated previous homography between the previous image and the first image of the two images. ,
Methods for stitching medical images.
According to claim 8,
The step of performing the relative evaluation,
A value obtained by subtracting the estimated horizontal shift value of the previous homography between the previous image and the first of the two images from the horizontal shift value of the estimated homography between the two images exceeds 0;
Homography is correct if the value obtained by subtracting the rotation value of the estimated previous homography between the previous image and the first of the two images from the rotation value of the estimated homography between the two images is less than the image rotation reference value. Judging by (1),
The value obtained by subtracting the estimated horizontal shift value of the previous homography between the previous image and the first of the two images from the horizontal shift value of the estimated homography between the two images is 0 or less, or the estimated horizontal shift value between the two images If the absolute value obtained by subtracting the rotation value of the previous homography estimated between the previous image and the first of the two images from the rotation value of the homography is equal to or greater than the image rotation reference value, the homography is set to not accurate (0). Including the step of judging,
Methods for stitching medical images.
According to claim 4,
The step of performing the final homography evaluation,
If the result of performing the self-evaluation or the result of performing the relative evaluation results in that the homography is correct (1) is output, it is determined that the homography estimated by the integrated feature detection algorithm is correct, and the result of performing the self-evaluation or the relative evaluation is correct. Including the step of determining that the homography estimated by the integrated feature detection algorithm is not correct when the result of performing the evaluation is all incorrect (0).
Methods for stitching medical images.
According to claim 1,
As a result of the evaluation, when it is determined that the homography estimated by the integrated feature detection algorithm is not accurate, re-estimating the homography between the selected two images using a local image-based estimation algorithm.
Methods for stitching medical images.
According to claim 11,
The step of re-estimating the homography,
selecting a corresponding first local area between the two images from each of the two images; and
Re-estimating a homography based on the selected first local area,
Methods for stitching medical images.
According to claim 12,
In the step of selecting the first local area,
Selecting a rectangular area corresponding to a length of a predetermined ratio of the size of each image as the first local area, centering on the position of the corresponding feature point between the two images in each of the two images.
Methods for stitching medical images.
According to claim 13,
The step of re-estimating the homography,
selecting a second local area included in a predetermined area of the two images from among the selected first local areas;
estimating a rotation component between two images based on the selected second local area;
estimating a motion component between the two images by applying the estimated rotation component; and
Re-estimating homography based on the estimated rotational component and movement component,
Methods for stitching medical images.
15. The method of claim 14,
The step of estimating the rotation component,
applying a Fourier transform to the selected local region; and
Estimating a rotation value by rotating the local area to which the Fourier transform is applied so that the local areas of the two images match.
Methods for stitching medical images.
According to claim 15,
The step of estimating the moving component,
correcting a second image of the two images through the estimated rotation component; and
Estimating a motion component between the two images based on a phase correlation with flatness algorithm,
Methods for stitching medical images.
17. The method of claim 16,
The step of re-estimating the homography,
estimating rotational components and performing correction reflecting rotational components for all corresponding local regions of the two images;
The flatness-based phase correlation with flatness algorithm is applied to all corresponding local regions of the two images for which the correction is completed, and one correspondence that has the largest maximum peak value and satisfies the preset flatness condition is obtained. deriving a local area to be; and
Comprising the step of re-estimating a homography by calculating a rotation component and a movement component based on the derived corresponding local area,
Methods for stitching medical images.
18. The method of claim 17,
The flatness,
Calculated as a smaller standard deviation value among the standard deviations of pixel values (intensities) of all corresponding pixels in the local area,
Methods for stitching medical images.
A computer-readable recording medium having recorded thereon a program that, when executed by at least one processor, causes the at least one processor to perform the method of any one of claims 1 to 18.
As a medical image stitching device,
a memory storing at least one instruction; and
Including at least one processor that executes the command in conjunction with the memory,
When the instruction is executed by the processor, it causes the processor to:
Select two consecutive images among the aligned images of the image set to be stitched,
Calculate a homography between the selected two images using an integrated feature detection algorithm based on feature points estimated based on a plurality of different feature detection algorithms;
Evaluate the accuracy of the calculated homography through a homographic evaluation algorithm,
Set to perform stitching between the selected two images by determining a homography between the two selected images according to the evaluation result.
Medical image stitching device.

Images