JP2013048696A - Image analysis device for ophthalmic disease, image analysis method for ophthalmic disease, and image analysis program for ophthalmic disease - Google Patents

Abstract

PROBLEM TO BE SOLVED: To perform automation and quantification of analysis so that the determination depending on the subjectivity of a reading doctor may be not caused as much as possible when performing diagnostic imaging.SOLUTION: An image analysis device for ophthalmic disease is characterized by comprising: an image input unit for acquiring three-dimensional voxel images of the head including the eyeballs; a voxel image extraction unit for extracting, from the input three-dimensional voxel images, at least one three-dimensional eyeball voxel image of the eye in an eyeball region separately for the right and left eyeballs; a projection image calculation unit for creating a first projection image obtained by projecting the extracted three-dimensional eyeball voxel image in an upward- and/or downward-direction; and a projection image analysis unit for setting an eye axis that links an ocular anterior pole apex and an ocular posterior pole apex on the first projection image, and calculating, as parameters for image analysis, at least the eye axis length between the ocular anterior pole apex and the ocular posterior pole apex, and the degree of left/right asymmetry on the basis of a region obtained by segmenting the posterior pole side to the right and left bordered by the eye axis.

Description

  The present invention relates to an image analysis technique for ophthalmic diseases, and more particularly to an image analysis technique for high myopia.

  Popular myopia in ophthalmological diseases is known as severe myopia, a severe illness that causes posterior polar tissue degeneration, eye axis elongation, strong myopia, and sometimes blindness. . Intensity myopia occurs due to degeneration or atrophy of the retina, the eyeball extends elliptically, and the axial length is 3.5 mm or more longer than normal vision. Of the approximately 67 million people over the age of 40 in the country, about 3.6 million people are considered to have high myopia. The retina is burdened, causing illnesses such as flying mosquito disease, retinal detachment, and glaucoma, narrowing the field of view, and may lead to vision loss and blindness.

  At present, the curative treatment for this disease is in the research stage such as Patent Document 1 and there is only symptomatic treatment that suppresses the symptoms and progression of the disease state by pharmacotherapy and surgical treatment, and early detection is important. However, since this disease gradually progresses over a long time span of 20 years or more, and in many cases it does not show an increase in axial length or a decrease in visual acuity at an early stage, current ophthalmologic diagnosis such as a fundus camera, etc. There was a problem that discovery tends to be delayed in the equipment.

  Myopia is determined on the basis of the value of the lens power D (diopter or diopter, D = 1 (m) / focal length (m)) required for refractive correction. In normal vision, the D value is 0, and in hyperopia, it is positive. A case where D is a negative value is myopia, and -6D or less is intense myopia. However, in the initial stage of intense myopia, the D value may not show an abnormal value, so that it is possible to diagnose earlier by using the characteristic that the eyeball extends in the front-rear direction.

  As a conventional diagnostic method for this disease, in addition to the above-described refraction inspection apparatus (Patent Document 2, “eye refractometer”, Canon, etc.) visual acuity measurement (Patent Document 3, “glasses, contact lens power determination system”, vision glasses) , Ocular axis measurement by ultrasound (Patent Document 4, “Ophthalmic Ultrasound Diagnostic Device”, Canon), Ocular axis measurement by optical interference (Patent Document 5, “Dimensional measurement device for living eye”, Topcon), Retina merge by fundus camera For example, there is an analysis of symptoms (Patent Document 6 “Analyzer for fundus disease”, Topcon) and an intraocular pressure measurement (Patent Document 7 “Non-contact tonometer”, Topcon) for diagnosis of complications with glaucoma.

Special table 2007-536272 gazette Japanese Patent Publication No. 5-054744 Japanese Patent No. 4014438 Japanese Patent Publication No. 3-058726 Japanese Patent No. 2723967 Japanese Patent No. 3479788 Japanese Examined Patent Publication No. 4-030294 JP 2007-121259 A

Brain J. Curtin, “The Posterior Staphyloma of Pathlogic Myopia”, Transaction of American Opthalmology Society, Vol. LXXXV, pp. 67-86, (1997).

  As for the pathological condition of intense myopia described in Non-Patent Document 1, 10 types classification based on a fundus camera image published by BJ Curtin in 1977 is known.

  In addition, it is possible to diagnose typical symptoms of “intensity myopia” such as extension of the eye axis and reduction of visual acuity using other conventional ophthalmologic diagnostic devices disclosed in Patent Documents 2 to 7, but degeneration of the posterior pole of the eyeball is possible.・ Regarding lesions such as hypertrophy, it was difficult to measure with conventional examination means approached from the front of the eyeball.

  In addition, when diagnosing degeneration / hypertrophy by introducing image diagnosis, there is a problem that qualitative determination that there is an asymmetry in either direction is likely to cause variation among interpretation doctors.

  The first object of the present invention is a stage before presenting typical symptoms of “intensity myopia” such as extension of the eye axis and decreased visual acuity, and degeneration / hypertrophy of the posterior pole of the eye, which is an early sign of the disease (referred to as posterior vineoma). Non-invasive and quantitative analysis. Then, the diagnosis of “intensity myopia” is confirmed early so that the progression of the disease state can be suppressed.

  A second object of the present invention is to realize automation and quantification of an analysis so that determination depending on the subjectivity of an interpreting doctor does not occur as much as possible in image diagnosis.

  The present invention introduces a popular three-dimensional tomography technique in the field of ophthalmology in the diagnosis of head disease, and acquires a three-dimensional image of the eyeball mainly by MRI (Magnetic Magnetic Resonance Imaging) of the head. We propose a method for noninvasive, quantitative and automatic analysis of the morphological feature parameters of the posterior pole of the eye, which correlates with the prognosis of this disease.

  For example, a three-dimensional MR voxel image of the left and right eyeballs cut out from the MRI image of the head using the technique described in Japanese Patent No. 3049021 (“Image Processing Device”, five majors) is disclosed in Japanese Patent No. 44799651 (“Medical Image Display Device”). ”, Hitachi Medical Corp.) and other techniques are used to prepare six-frame projected images that have been volume-rendered in six directions, up, down, left and right, and to analyze the eyeball.

  Subsequently, an eye axis is set for each of the four projected images in the vertical and horizontal directions, and a geometrical symmetry in the horizontal and vertical directions in the radial area from the center of gravity to the posterior pole, The sharpness of the rear curved surface is calculated. Furthermore, as an option, an analysis is performed as to whether the number of rear projections (vines) is single or plural from the rear projection image.

  The setting of the eye axis is based on automatic extraction, and in some cases, correction is made by interactive means. All other calculation processes are automated so that the parameters calculated among operators do not vary. Based on these 4 + option 1 items of axial length, left-right symmetry, vertical symmetry, sharpness, and number of protrusions, we classify the pathological conditions of high myopia and use them for prognosis determination and treatment policy formulation.

  According to one aspect of the present invention, an image input unit that acquires a three-dimensional voxel image of a head including an eyeball, and 3 of the eye region of at least one of eyeball regions by left and right eyeballs from the input three-dimensional voxel image. A voxel image extraction unit that extracts a three-dimensional eyeball voxel image; a projection image calculation unit that creates a first projection image obtained by projecting the extracted three-dimensional eyeball voxel image in at least one of the upper and lower directions; and the first projection. Set an eye axis connecting the anterior ocular vertex and the posterior ocular vertex with respect to the image, and at least the ocular axial length between the ocular anterior vertex and the ocular posterior vertex, with the eye axis as the boundary An image analysis apparatus for ophthalmic diseases, comprising: a projection image analysis unit that calculates a degree of asymmetry in the left-right direction based on a region obtained by dividing the posterior pole side into left and right, and a parameter for image analysis Provided.

  According to the present invention, intensity myopia can be determined based on the length of the axial axis or the degree of asymmetry in the left-right direction. The parameter includes a parameter for specifying the intensity myopia.

The three-dimensional voxel image is acquired using MRI.
The projection image calculation unit creates a second projection image obtained by projecting the three-dimensional eyeball voxel image in at least one of the left and right directions, and the image analysis unit generates an eye axis with respect to the second projection image. And at least the degree of asymmetry in the vertical direction based on the axial length between the apex of the front of the eyeball and the apex of the back of the eyeball and the area obtained by dividing the back pole side up and down with the eye axis as the boundary. , And is calculated as a parameter for image analysis. Thereby, the asymmetry in the vertical direction can be calculated.

  In addition, the image analysis unit calculates a sharpness degree of the arcuate curvature of the back pole as an image analysis parameter for at least one of the first and second projection images. . Thereby, the back pole sharpness can be calculated.

  In addition, the image analysis unit calculates a value obtained by averaging the parameters of at least one of the left, right, and upper and lower sets for at least one of the first and second projection images. To do. Thereby, the average value of asymmetry and rear pole sharpness can be calculated.

  Further, the projection image calculation unit creates a third projection image projected backward with respect to the three-dimensional eyeball voxel image, and the image analysis unit performs posterior to the third projection image. The number of peak points corresponding to the pole protrusions is calculated as a parameter for image analysis. Thereby, the number of rear pole protrusions can be calculated.

  In addition, the image analysis unit extracts an outline of an eyeball part from the third projection image, sets a center of gravity of the eyeball part, and performs image analysis on the calculated position of each peak point. As a parameter, the distance from the center of gravity, which is a polar coordinate based on the center of gravity point, and the angle from the center of gravity to the peak point direction are calculated. Thereby, the position of the number of rear pole protrusions can be calculated.

  Further, the image analysis unit extracts an outline of an eyeball from the projection image, the anterior ocular vertex and the ocular posterior apex are located on the outline of the eyeball, and the eyeball The straight line connecting the front pole apex and the posterior pole apex of the eyeball is equal in distance to the intersection with the contour line extending vertically from left and right or up and down from several vertices of the straight line, The eye axis is automatically set by setting the apex of the front of the eyeball and the apex of the back of the eyeball. Thereby, it is possible to automatically set the eye axis.

  Alternatively, the image analysis unit extracts an outline of an eyeball part from the projection image, further sets a center of gravity of the eyeball part, and the eyeball front pole vertex and the eyeball rear pole vertex are the eyeball part. And the posterior pole of the eyeball and the posterior part of the eyeball so that a polygonal line connecting the three points of the apex of the anterior eyeball, the center of gravity and the posterior pole of the eyeball is as straight as possible. The eye axis may be automatically set by setting the pole apex and the barycentric point.

  Further, the image analysis unit extracts an outline of an eyeball part from the first or second projection image, and displays a screen for either the frontal vertex of the eyeball or the vertex of the back of the eyeball. A straight line connecting the anterior ocular pole vertex and the posterior pole apex with respect to the anterior ocular pole vertex or the posterior pole apex set by the interactive apex input unit set in an interactive form with A vertex correction unit that automatically corrects the intersection with the contour line on the extended line is provided. Thereby, it is possible to manually set the eye axis.

Further, the image analysis unit extracts an outline of the eyeball part from the first or second projection image, further sets a measurement center point on the eyeball part, and sets the posterior pole of the eyeball from the measurement center point. A straight line A and a straight line B in the direction of the first angle and the first angle are set around the center of gravity for measurement with respect to the reference line connecting the vertices, and are surrounded by the straight line A and the reference line. Further, an area A inside the contour line and an area B inside the contour line surrounded by the straight line B and the reference line are calculated, and the ratio of the area A to the area B is given as the degree of asymmetry. It is characterized by that. As a result, asymmetry can be calculated.
The first angle is, for example, not less than 40 degrees and less than 90 degrees.

  The image analysis unit extracts an outline of an eyeball from the first or second projection image, further sets a measurement center of gravity in the eyeball, and determines the posterior pole apex of the eyeball from the measurement center of gravity. A straight line A and a straight line B are set in the direction of the second angle and the second angle around the center of gravity for measurement with respect to the reference line connecting the straight line A and the intersection A between the straight line A and the contour line An angle formed by three vertices of the posterior pole apex of the eyeball and the intersection B of the straight line B and the contour line is given as the sharpness degree. Thereby, the back pole sharpness can be calculated. Note that the second angle is, for example, about 10 degrees to 40 degrees.

  Alternatively, the image analysis unit may define a straight line C in the direction of the + second angle and the −second angle with respect to the reference line connecting the vertex of the posterior pole of the eyeball from the center of gravity for measurement with respect to the center of gravity for measurement. A straight line D is additionally set, and an area C inside the contour line surrounded by the straight line C and the reference line, and an area D inside the contour line surrounded by the straight line D and the reference line are added. It may be calculated and a larger value of the ratio between the area C and the area D and the ratio between the area A and the area B may be set as the degree of asymmetry.

  In addition, the image analysis unit, as the center of gravity for measurement, on the eye axis that connects the apex of the front of the eyeball and the apex of the back of the eyeball, a point separated from the apex of the back of the eyeball by a predetermined distance It is characterized by giving. Thereby, the center-of-gravity point for measurement can be set.

  For example, the image analysis unit may set 1/2 (for example, 87 pixels, 13.4 mm) of the average eye axis length of a normal eyeball as the predetermined distance.

  According to another aspect of the present invention, an image input step for acquiring a three-dimensional voxel image of a head including an eyeball, and at least one eye of an eyeball region for each of right and left eyeballs from the input three-dimensional voxel image. A voxel image extraction step of extracting a three-dimensional eyeball voxel image, a projection image calculation step of creating a first projection image obtained by projecting the extracted three-dimensional eyeball voxel image in at least one of the upper and lower directions, and the first Set the eye axis connecting the anterior-polarity vertex and the posterior-polarity vertex to the projected image, and at least the eye-axis length between the anterior-polarity vertex and the posterior-polarity vertex, and the eye axis And a projected image analysis step for calculating a degree of asymmetry in the left-right direction based on a region obtained by dividing the posterior pole side into left and right sides, and an image of the ophthalmic disease characterized by comprising:析方 there is provided a method.

  The present invention may be a program for causing a computer to execute the above-described ophthalmic disease image analysis method, or may be a computer-readable storage medium for recording the program.

  The present invention proposes a novel pathological classification method based on shape feature parameters that are highly correlated with clinical prognosis by non-invasively measuring three-dimensional eye features of intense myopia using MRI. To do. Based on the 3D MR projection image rendered from 6 directions, the axial length, the left / right and up / down symmetry of the back of the eyeball, and the sharpness of the back of the eyeball, which are likely to be highly correlated with complications, are automatically calculated. Can be classified based on these conditions.

  According to the present invention, it is possible to contribute to an early diagnosis of “intensity myopia” without imposing a burden on the patient, and to provide objective diagnosis data that does not depend on the subjectivity of the interpretation doctor.

  Further, image processing can be realized with inexpensive general-purpose personal computer software, and a special graphics board such as a GPU (Graphics Processing Unit) is not required, so that an inexpensive diagnostic service can be provided.

It is a flowchart figure which shows the flow of the rough diagnostic process of an intensity | strength myopia diagnosis. It is a flowchart figure which shows the outline | summary of the intensity | strength myopia diagnostic method using the image analysis technique of the ophthalmic disease using head MRI imaging | photography data by this Embodiment. 2B (a) shows an example of an MRI projection image taken from above, FIG. 2B (b) is an eyeball projection image viewed from the left, FIG. 2B (c) is an eyeball projection image viewed from above, and FIG. 2B (d) is a projected image of the eyeball viewed from behind. FIG. 3A is a functional block diagram illustrating a configuration example of an image analysis apparatus according to an embodiment of the present invention, and FIG. 3B is a functional block diagram illustrating a configuration example of a projection image analysis unit. is there. It is a figure which shows the measurement group of intensity myopia using the image acquired by MRI by this Embodiment, for example, it calculates | requires from the image of (b) to (d) of FIG. It is. It is a flowchart figure which shows an example of the intensity | strength myopia diagnostic method by MRI by this Embodiment. These are front and rear, right and left images (300 × 300 pixels), and the front and rear projection images are not essential for analysis, but are acquired for reference. It is a figure which shows the detail of the measurement parameter of the intensity | strength myopia by MRI. It is a figure which shows the detail of the measurement parameter of the intensity | strength myopia by MRI. It is a figure which shows the analysis process of the projection image of an up-down direction or the left-right direction of step S31-1. It is a figure which shows the analysis process of the projection image of an up-down direction or the left-right direction of step S31-1. It is a figure which shows the analysis process of the projection image of an up-down direction or the left-right direction of step S31-1. It is a figure which shows the mode of the back projection image process from step S51 of FIG. FIG. 12A is a diagram illustrating processing subsequent to FIG. FIGS. 12B to 12E are diagrams showing the processing. It is a figure explaining the detail of the binarization process of an image, and an outline extraction process (S32). FIG. 14A is a flowchart showing the flow of the neighborhood region integration processing, and FIGS. 14B to 14D are views showing the processing. It is a figure which shows the outline | summary of the algorithm of the automatic setting of an eye axis. It is a figure which shows the detail of the automatic setting process (step S33) of an eye axis. It is a figure which shows the flow of the process following FIG. It is a figure which shows the correction | amendment of a back pole side edge point, the determination process of an eye axis, and the mode. It is a figure which shows the flow of the correction process of an eye axis, and a figure which shows the mode. It is a figure which shows the detail of the creation process (S38) of the edge point in a back pole outline part, and its mode. FIG. 20A is a flowchart showing the flow of the smoothing filter process (FIG. 11A: S52). FIG. 20B is a diagram illustrating an example of the smoothing filter. It is a flowchart figure which shows the flow of the extraction process (refer FIG.11 (b): S53) of the light / dark peak point with respect to a projection image. It is a flowchart figure which shows the flow of the clear process of the surrounding pixel of the peak point extracted by the process of FIG. 21 (step S54 of FIG. 11). It is a flowchart figure which shows the calculation process of the gravity center point with respect to an original image (FIG. 12 (a): step S55). It is a flowchart figure which shows the detail of the binarization process (step S32) of a pixel. FIG. 25A is a flowchart showing the flow of neighboring pixel connection processing and region extraction processing (see step S32-3, FIG. 13D). FIG. 25B is a diagram illustrating four of the eight neighboring pixels. FIG. 26A is a flowchart showing the flow of the neighborhood region integration process (see FIG. 14, step S32-4, FIG. 14B). FIG. 26B is a diagram illustrating an example of 8 neighboring pixels. It is a figure which shows the flow of the process following FIG. FIG. 15 is a flowchart showing a minute area deletion process (step S32-5) for deleting a minute area such as dust around and reducing an error in determining an eye axis (see FIG. 14C). ). It is a flowchart figure which shows the flow of an outline extraction process (refer FIG.14 (d)). It is a flowchart figure which shows the flow of a calculation process of a front pole side and a back pole side edge point, and a gravity center point (S33-1: FIG. 15B). It is a flowchart figure which shows the correction | amendment and determination process of a front pole side edge point (refer S33-3 and FIG. 16). It is a flowchart figure which shows the flow of the distance difference calculation process of the front-and-back edge point and gravity center direction which is a process of step S33-3-1 of FIG. 31 (FIG. 16, step S33-4). It is a flowchart figure which shows the flow of correction | amendment and a determination process (S33-4) of a gravity center point (refer FIG.16 (c)). It is a flowchart figure which shows the correction | amendment of a back pole side edge point, and the determination process of an eye axis (automatic axis setting process S33 of FIG. 16). It is a flowchart figure which shows the flow of the determination process of the linearity of 3 vertices (refer FIG. 18). It is a flowchart figure which shows the automatic correction | amendment process (S35-3) to the outline of a manually set edge point. It is a flowchart figure which shows the flow of the process (refer step S38-1, Fig.19 (a)) which adds a 90 degree and -90 degree edge point to the direction orthogonal to the already determined eye axis. It is a flowchart showing the flow of additional processing of edge points of theta _w and - [theta] _w (± first angle) (see step S38-2, FIG 19 (a), (c)). It is a flowchart figure which shows the flow of an addition process of the edge point (+/- 2nd angle) of (theta) _n and-(theta) _n degree (refer step S38-3, FIG. 19 (a), (d)). (A) is a figure which shows the calculation method of the circular arc area in the case of a left direction projection image, (b) is a figure which shows the ratio of the area calculated | required by (a). FIG. 41 is a flowchart illustrating an example of calculation processing of an arc area in a rear pole portion described in FIG. 40. (A) is a figure which shows the calculation method (narrow area | region) of the circular arc area in a rear pole part, (b) is a figure which shows an area ratio. FIG. 43 is a flowchart illustrating an example of calculation processing of an arc area in a rear pole portion described in FIG. 42. It is a figure explaining the calculation method of angle (eta) which the back pole part 3 edge point (sharpness) comprises. It is an example of a six-direction projection image to be analyzed, and it can be seen that projection images from six directions, i.e., left and right, up and down, and front and rear, are obtained. It is a figure which shows the example of the setting screen of the analysis parameter in analysis. It is a figure which shows an example of an analysis result, and is a figure which sees the analysis result based on an analysis process on one screen. It is a figure corresponding to FIG. 45, and is a figure which shows an example of the confirmation screen of the auxiliary line setting condition at the time of analysis. It is a figure which shows the example of the 6 direction projection image of the analysis object in the case of intense myopia. It is a figure which shows the example of a setting screen of an analysis parameter. It is a figure which shows the example of the result of having analyzed in the same setting as FIG. It is a figure which shows the confirmation screen of the auxiliary line setting condition at the time of analysis. It is a figure which shows the example of the 6-direction projection image of the 6-direction projection image (mild) of the analysis object. It is a setting screen of the analysis parameter which made the left eye an analysis symmetry. It is a figure which shows the analysis result based on FIG. It is a figure which shows the example of the confirmation screen of the auxiliary line setting condition at the time of analysis. It is a confirmation screen of the auxiliary line setting situation at the time of analysis, and is a diagram showing a state in which correction of the eye axis (manual correction in this example) is added. This is the analysis parameter setting screen. For example, the setting value is not changed from the setting screen shown in FIG. 54, but the setting value can be confirmed on the screen. It is a figure which shows an example of the reanalysis result after correction.

  Intensity myopia inspection technology according to an embodiment of the present invention uses three-dimensional tomography technology to acquire projection images from a plurality of viewpoints including an eyeball, and performs diagnosis of intensity myopia based on the projection images. It is characterized by.

  Hereinafter, an image analysis technique for an ophthalmic disease using head MRI imaging data according to an embodiment of the present invention will be described in detail with an example of intensity myopia.

  First, the flow of a rough diagnosis process for intensity myopia diagnosis will be described with reference to the flowchart of FIG. First, in step S1, diagnosis is performed using a general ophthalmologic diagnosis device. Common ophthalmic diagnoses include visual acuity tests, refractive index tests, retinal lesion tests, visual field tests, and the like. In step S2, it is determined whether or not the diagnosis is confirmed. If the diagnosis is confirmed, the treatment policy is also confirmed (to step S8).

  In the case of No in step S2, head MRI imaging and eyeball part rendering (separate left and right eyeballs) are performed. In step S4, projection image data of the head in six directions, up, down, left and right, front and rear, is generated for each left and right eyeball. Next, image analysis processing is performed on a two-dimensional (2D) projection image. In step S6, shape parameter and intensity myopia type classification processing is performed. Based on the result of step S6, the treatment policy is examined (step S7), and the process proceeds to step S8 to determine the treatment policy.

  FIG. 2A is a flowchart showing an overview of an intensity myopia diagnosis method using an image analysis technique for ophthalmic diseases using head MRI imaging data according to the present embodiment. FIG. 2B (a) shows an example of a tomographic image obtained by slicing an eyeball part from an MRI obtained by imaging the head of a human body, and for a three-dimensional voxel image obtained by extracting one eyeball part from the head MRI imaging data. 2B (b) is a projected image of the eyeball viewed from the left, FIG. 2B (c) is a projected image of the eyeball viewed from above, and FIG. 2B (d) is a projected image of the eyeball viewed from the rear.

  As shown in FIG. 2A, first, in step S11, head MRI imaging is performed, and voxel data of the eyeball region is extracted for each of the left and right eyeballs (the method described in Japanese Patent No. 3049021 can be used).

  Next, in step S21, a two-dimensional projection image in the vertical direction or the horizontal direction is volume-rendered (one or two in the vertical direction or the horizontal direction) for each eyeball voxel data separately on the left and right sides to create a projection image. (The method described in Japanese Patent No. 4497965 can be used).

  Next, in step S31, an eye axis is set for each of 2 to 4 projection images in the vertical direction or the horizontal direction for the left and right eyeballs. Optionally, the shape characteristic of the back pole is measured from the projection image.

  Further, in step S51, a rearward two-dimensional projection image is rendered for each eyeball voxel data for each of the left and right eyeballs to create a projection image (one each). Using this projection image, the number of protrusions is arbitrarily measured.

  FIG. 2B (a) is a tomographic image obtained by a head MRI of a healthy person. The projection image of the eyeball seen from the left of FIG. 2B (b), the projection image of the eyeball seen from the upper part of FIG. 2B (c), and the projection image of the eyeball seen from the rear of FIG. 2B (d) are used.

  FIG. 3A is a functional block diagram illustrating a configuration example of an image analysis apparatus according to an embodiment of the present invention, and FIG. 3B is a functional block diagram illustrating a configuration example of a projection image analysis unit. is there.

  The image analysis apparatus shown in FIG. 3A includes an image input unit 1 that inputs an image acquired by MRI, a voxel image extraction unit 3 that extracts voxel data from the image, and a projection image that calculates a projection image from the voxel data. The calculation unit 5 includes a projection image analysis unit 7 that performs image analysis of the projection image, and an output unit 11 that outputs an analysis result. As shown in FIG. 3B, the image analysis unit 7 includes an asymmetry calculation unit 7-1 that calculates asymmetry, a sharpness degree calculation unit 7-2 that calculates the degree of sharpness, and a peak of the rear pole protrusion. A rear pole protrusion peak calculation unit 7-3 for calculating the eye axis, an eye axis automatic setting unit 7-4 for automatically setting the eye axis, and a vertex correction unit 7-5 for correcting the vertex.

FIG. 4 is a diagram showing intensity myopia measurement parameters using projection images acquired by MRI according to the present embodiment, and is obtained from, for example, the projection images from (b) to (d) of FIG. 2B. FIG. 4A is a diagram showing a method of obtaining the axial length L1 that is the length of the eye along the axial axis from the vertical and horizontal projection images. For example, the axial length L1 is calculated in the vertical and horizontal directions. The average value of the measurement results from the four projected images is determined. In the projection image of FIG. 2B (b), the horizontal direction from the left end portion of the small semicircle region (crystal) on the left side to the right end portion of the eyeball in the eyeball portion displayed in white (FIG. 2B (b) In the case of, find the length of the line stretched slightly upwards). FIG. 4B is a diagram showing a method for obtaining the posterior pole sharpness θ, and two intersections with the posterior polar contour line that extends from the measurement center of gravity to the posterior pole direction at an angle of ± θn with C1 as the eye axis. The angle θ (rear pole sharpness θ) formed by the lines L ₂ and L ₃ connecting the vertebra and the apex of the back of the eyeball is determined by the average value of the measurement results of the four projected images in the vertical and horizontal directions. FIG. 4C is a diagram illustrating a method for obtaining the left-right symmetry. For example, the left-right symmetry is determined by the average value of the measurement results obtained from the two vertical projection images. Specifically, in FIG. 4 (c), the area E and the nose of the posterior pole divided by the eye axis C1 (the right eye: 0 to + θ _W , the left eye: 0 to −θ _W ). The ratio with the area N of the rear pole of the side region (for the right eye: 0 to −θ _W , for the left eye: 0 to + θ _W ) is obtained. Actually, as will be described later, the ratio of the area E to the area N is similarly obtained using θ _n smaller than θ _W , and the one having the larger asymmetry degree of the area ratio is adopted. This symmetry is an important parameter having a high correlation with the severity. FIG. 4D is a diagram illustrating a method of obtaining vertical symmetry. For example, the upper region T (0 to + θ _W ) and the lower region B (0 to −) whose vertical symmetry is partitioned by the eye axis C1. The ratio T / B area ratio with θ _W ) is obtained, and the average value of the upper and lower values is obtained. Actually, as will be described later, the ratio between the area T and the area B is similarly obtained using θ _n smaller than θ _W , and the one with the larger area ratio is adopted. FIG. 4 (e) obtains the number of rear pole protrusions from the projected image in the rear direction, and obtains whether or not there are a plurality of rear vines.

Table 1 shows an example of type classification of intensity myopia by MRI. As shown in Table 1, in the four-digit type classification code (all 54 patterns from 0000 to 2221), the first place is the axial length, normal is “0”, and abnormality is “1”. Position 10 is a horizontal direction of symmetry (determined at the posterior pole lateral direction theta _w degree arc area ratio above and below the projection image), 0: symmetry, 1: Nose closer (area ratio less than 90%), 2: Close to the ear (area ratio of 110% or more). 100's place are vertical symmetry, it determines at best under the direction theta _w degree arc area ratio rear left and right images. 0: Symmetric, 1: Lower (less than 90%), 2: Upper (110% or more). The position of 1000 is the rear pole sharp type, which is determined by the angle formed by the center of the rear pole tip in the upper, lower, left and right projected images. 0: Intermediate, 1: Spindle type (less than 140 degrees), 2: Barrel type (150 degrees or more). With this four-digit code, it is possible to classify and organize the pathological conditions related to intense myopia.

  FIG. 5 is a diagram showing an example of an intensity myopia diagnostic technique based on MRI according to the present embodiment. FIG. 6 shows projection images (300 × 300 pixels) in the front-rear, left-right, up-down directions of the eyeball. The projection images in the front-rear direction are not essential for analysis, but are acquired for reference.

  First, in step S11, head MRI imaging is performed (voxel data of the eyeball region is extracted for each of the left and right eyeballs). As the processing algorithm, for example, the method described in Japanese Patent No. 3049021 can be used. For each eyeball voxel data for each left and right eyeball, two-dimensional projection images in six directions, up and down, left and right, and front and rear, are rendered (total of six images). As the processing algorithm, for example, the method described in Japanese Patent No. 4497965 can be used.

  More specifically, for the four projected images in the upward, downward, leftward, and rightward directions, the eye axis is set automatically or manually by a method described later, and parameters are calculated. The rear projection image is for automatically or visually determining the number of protrusions of the posterior vine, and when there are a plurality of projections, the result of left-right symmetry or vertical symmetry described later is not properly determined. It may be used to supplement the results. As a method for automatically extracting the number of protrusions of the posterior vine, it can be realized by an algorithm for automatically extracting peak points from the nuclear medicine image described in Patent Document 8. Although the forward projection image is not subject to analysis, it is used as a reference for the interpretation doctor to determine whether or not the automatically set eye axis is appropriate.

  Subsequently, as shown in FIG. 7A (b), which will be described later, a measurement gravity center Po is defined at a predetermined distance from the rear pole edge point Pb on the eye axis set for each of the four projected images. The center of gravity Po for measurement is close to the midpoint of the eye axis Pf-Pb or the center of gravity of the eyeball Pg in the case of a normal eye, but is biased toward the rear side of the eyeball in high myopia. Based on these set eye axes and measurement center of gravity, the four parameters shown in FIG. 4 are analyzed as follows.

Next, in step S31-1, the eye axis is set for the four projected images in the vertical and horizontal directions, and the shape characteristic of the posterior pole is measured (repeated by four projected images, but at least one projection is performed). Processing is possible if there is an image). Next, in Step S31-2, the left-right direction symmetry, the up-down direction symmetry, and the rear pole sharpness are calculated based on the measurement results of the four projected images. Here, the symmetry in the horizontal direction calculates the average value of the measurement results of the vertical projection image, the symmetry in the vertical direction calculates the average value of the measurement result of the horizontal projection image, and the rear polar sharpness is An average value of the measurement results of the four projected images is calculated. Next, in step S51, the number of protrusions is measured for the backward projected image.
With the above processing, the parameters shown in FIG. 4 can be obtained for the actual eyeball.

  7A and 7B are diagrams illustrating specific examples of the morphological analysis parameters for intensity myopia. Based on FIGS. 7A and 7B, an outline of an analysis method of each parameter will be described.

(1) Axial length: The axial length is described based on FIG. 7A (a). The lengths of the eye axes P _f -P _b set for each of the projected images in the four directions of the upward direction, the downward direction, the left direction, and the right direction are calculated individually. Subsequently, the angle α formed between the eye axis P _f −P _{b of} the left or right projected image and the X axis (horizontal line) is calculated, and the eye axis P _f −P of the upward or downward projected image. _An angle β formed by _b and the Y axis (vertical line) is calculated. For the axial length calculated from the projected image in the left or right direction, the axial length calculated from the projected image in the upward or downward direction is obtained by multiplying 1 / cosβ to correct the tilt in the vertical direction. Is multiplied by 1 / cos α to correct the horizontal inclination. However, since the correction amount by these correction coefficients is very small, there is no effect even if omitted in practice. Finally, an average value ((L ₁ + L ₂ + L ₃ + L ₄ ) / 4) of the axial lengths L ₁ , L ₂ , L ₃ , and L ₄ calculated based on the projection images in the four directions is a predetermined value. It is determined as 0 (normal) or 1 (abnormal) depending on whether a threshold value (eg, 175 pixels) is exceeded.

により算出できる。そして、４方向の画像を基に算出された後極尖鋭度θ_１、θ_２、θ_３、θ_４の平均値（（θ_１＋θ_２＋θ_３＋θ_４）／４）が、２つの下限しきい値（例．140度）と上限しきい値（例．150度）を設定し、両しきい値範囲であれば０（中間型）、下限しきい値未満は１（紡錘型）、上限しきい値を超えたら２（樽型）と判定する。 (2) Rear pole sharpness: The rear pole sharpness will be described based on FIG. 7A (b). Calculating a curvature at Kokyoku edge point P _b every four directions of each projected image upward, downward, leftward and right direction. However, since a difference in solids hardly occurs at the radius of curvature, the rear pole of a straight line extending in the direction of a predetermined angle (for example, ± θ _n degrees) from the center of gravity P ₀ to the left or right or up and down around the rear pole edge point P _b. the intersection of the contour and _{P + n} and _{P -n,} calculates the angle _P + _{n -P} b -P _-n three vertices of the two intersections and Pb, the angle (360 ° increments) the posterior pole acuity It is defined as In particular,

Can be calculated. Then, the average value ((θ ₁ + θ ₂ + θ ₃ + θ ₄ ) / 4) of the rear polar sharpness θ ₁ , θ ₂ , θ ₃ , θ ₄ calculated based on the images in the four directions has two lower limits. Set threshold value (eg 140 degrees) and upper threshold value (eg 150 degrees), 0 (intermediate type) if both threshold ranges, 1 (spindle type) if lower threshold value, upper limit If it exceeds the threshold, it is determined to be 2 (barrel type).

続いて、右眼の場合はＳ_＋ｗ×100／Ｓ_−ｗを広領域Ｅ／Ｎ面積比率およびＳ_＋ｎ×100／Ｓ_−ｎを狭領域Ｅ／Ｎ面積比率とし、左眼の場合はＳ_−ｗ×100／Ｓ_＋ｗを広領域E/N面積比率およびＳ_−ｎ×100／Ｓ_＋ｎを狭領域E/N面積比率として算出する。そして、２種類の広領域E/N面積比率と狭領域E/N面積比率より100との差が大きい方のE/N面積比率を左右対称性の判定に採用し、上方向と下方向の投影画像を基に算出された２投影画像のE/N面積比率（耳側／鼻側比率）［％］の平均値を左右対称性とする。最後に、２つの下限しきい値（例．９０％）と上限しきい値（例．１１０％）を設定し、両しきい値範囲であれば０（対称）、下限しきい値未満は１（鼻側非対称）、上限しきい値を超えたら２（耳側非対称）と判定する。 (3) Left-right symmetry: Right-left symmetry will be described with reference to FIG. 7B (c). Straight lines extending in the direction of two types of angles (for example, ± θ _w degrees and ± θ _n degrees) from the measurement center of gravity P _{0 of the} upward or downward projection image to the left and right around the rear pole edge point P _{b Let} P _{+ w} and P _−w and P _{+ n} and P _−n be the intersections on the rear polar outline. The area of the arc surrounded by P ₀ -P _b -P _{+ w} is S _{+ w} , the area of the arc surrounded by P ₀ -P _b -P _-w is S _-w , and P ₀ -P _b -P _{+ n.} The arc area of four regions is calculated with the area of the arc surrounded by S _{+ n} and the area of the arc surrounded by P ₀ −P _b −P _−n being S _−n .

Subsequently, in the case of the right eye, S _{+ w} × 100 / S _−w is the wide region E / N area ratio and S _{+ n} × 100 / S _−n is the narrow region E / N area ratio, and in the case of the left eye, S _{− w} × 100 / S _{+ w} is calculated as the wide region E / N area ratio, and S _−n × 100 / S _{+ n} is calculated as the narrow region E / N area ratio. Then, the E / N area ratio, which has a larger difference between the two types of wide area E / N area ratio and 100 than the narrow area E / N area ratio, is adopted for the determination of left-right symmetry. The average value of the E / N area ratio (ear side / nose side ratio) [%] of the two projected images calculated based on the projected images is defined as left-right symmetry. Finally, two lower threshold values (eg 90%) and upper threshold value (eg 110%) are set. If both threshold ranges are 0 (symmetric), less than the lower threshold is 1 (Nose side asymmetry), 2 (ear side asymmetry) is determined if the upper threshold is exceeded.

  The reason for calculating the E / N area ratio in the two large and small angle ranges is that the rear pole edge point Pb is slightly left and right (usually the ear side) when the expansion volume of the rear pole is relatively small. This is due to the characteristic of the myopia progression of intense myopia that the rear pole part is shifted to the left and right (usually the nasal side) when the rear volume part is displaced and the expansion volume of the rear pole part is relatively large. In the former case, the narrow area E / N area ratio is determined to be asymmetric, but the wide area E / N area ratio is easily determined to be symmetric. On the other hand, in the latter case, the narrow area E / N area ratio is determined to be asymmetric (usually the ear side), and the wide area E / N area ratio is also determined to be asymmetric in the opposite direction (usually the nose side). The latter asymmetry was found to be more prominent.

(4) Vertical symmetry: Vertical symmetry will be described with reference to FIG. 7B (d). From the center of gravity Po for each measurement of the left or right projection image, two angles (for example, ± θ _w degrees and up and down) with respect to the back pole edge point Pb up and down are used in the same manner as the above-described left-right symmetry calculation method. Intersections on the rear polar outline of the straight line extending in the direction of ± θ _n degrees) are defined as P _{+ w} and P _−w and P _{+ n} and P _−n . The area of the arc surrounded by P ₀ -P _b -P _{+ w} is S _{+ w} , the area of the arc surrounded by P ₀ -P _b -P _-w is S _-w , and P ₀ -P _b -P _{+ n.} The arc area of the four regions is calculated with the area of the arc surrounded by S _{+ n} and the area of the arc surrounded by P ₀ −P _b −P _−n being S _−n . Subsequently, S _{+ w} × 100 / S _−w is calculated as a wide region T / B area ratio, and S _{+ n} × 100 / S _−n is calculated as a narrow region T / B area ratio. Then, the T / B area ratio having the larger difference between the two types of wide area T / B area ratio and 100 than the narrow area T / B area ratio is adopted for the determination of the vertical symmetry. The average value of the T / B area ratio (upper / lower ratio) [%] of the two projection images calculated based on the projection images is defined as vertical symmetry. Finally, two lower threshold values (eg 90%) and upper threshold value (eg 110%) are set. If both threshold ranges are 0 (symmetric), less than the lower threshold is 1 (Lower asymmetric), if the upper threshold value is exceeded, 2 (upper asymmetric) is determined. However, in the case of high myopia, it is rare to be determined as upper asymmetric.

  When the above analysis is performed, the axial length is 0 (normal) or 1 (abnormal), and the posterior sharpness is 0 (intermediate type), 1 (spindle type), 2 (barrel type), Left / right symmetry is 0 (symmetric), 1 (nasal asymmetry) or 2 (ear asymmetry), and vertical symmetry is 0 (symmetry), 1 (lower asymmetry), 2 (upper asymmetry). The pathological conditions of high myopia are theoretically classified into 54 types including normal eyes.

  Below, the content of the specific pre-processing with respect to a projection image for calculating the form analysis parameter mentioned above is demonstrated.

  FIG. 8A, FIG. 9A, and FIG. 10A are diagrams illustrating the analysis processing of the projected image in the vertical direction or the horizontal direction in step S31-1. FIG. 8B, FIG. 9C, and FIG. 10D are diagrams showing an outline of an image corresponding to each process.

First, in step S32, the binarization process and contour extraction of the projected image are performed on the projected image in either the horizontal direction or the vertical direction (see FIG. 8B). Next, in step S33, an automatic setting process for the eye axis C is performed (step S33). In Step S34, it is determined whether or not the set eye axis is appropriate. If Yes, the process proceeds to Step S36, and if No, the process proceeds to Step S35. In step S35, an eye axis correction process is performed (see FIG. 8D). Here, the intersection of the eye axis C and the outline of the eyeball is moved up and down to correct the eye axis to C ′.
As shown in FIG. 9 (a), in step S36, the axial length (distance Da from the front pole side edge point to the rear pole side edge point) is measured with the determined eye axis C (FIG. 9 (b)). ^{: Da = {(X2-X1} ) 2 + (Y2-Y1) 2} 1/2).

  In step S37, a measurement center of gravity Po (X0, Y0) is added to the eye axis C at a predetermined distance position (normal eye axis length / 2) from the rear pole side edge point. Here, X0 = (X1-X2) Ds / Da + X2 and Y0 = (Y1-Y2) Ds / Da + Y2 (see FIG. 9C). In the embodiment, Ds = 87 pixels (13.4 mm).

In step S38, edge points P _{+ w} , P _{+ n} , P _−w , and P _−n in the direction of positive / negative θ _w · θ _n degrees from the center of gravity Po for measurement centering on the rear pole side edge point are formed on the rear pole side contour. It adds (refer FIG.9 (d)).

Next, as shown in FIG. 10, in step S39, the arc area ratio is calculated in the range of the rear pole portion -θ _n (for example, 22.5 degrees is set) to 0 degrees and 0 degrees to + θ _n degrees (narrow region symmetry). (See FIG. 10B). In step S40, the calculation of the posterior pole portions - [theta] _w (e.g. setting the 45 °) to 0 ° and 0 ° ~ + theta _w arc area ratio in the range (wide area symmetry) to (see FIG. 10 (c)). In step S41, the larger one of the narrow area area ratio and the wide area area ratio is set as the symmetry parameter. Next, in step S42, an angle formed by connecting the three edge points of −θ _n degrees, 0 degrees, and + θ _n degrees of the rear pole contour is calculated (see FIG. 10D: sharpness θ).

  FIG. 11A is a flowchart showing a flow of backward projection image processing from step S51 of FIG. First, in step S52, smoothing filter processing (noise removal) is performed on the projected images in the backward direction (one each on the left and right) (see FIG. 11B). Next, in step S53, one shade peak point Pn (n = 1, 2,...) Is first extracted from the updated projection image. FIG. 11C is a diagram illustrating a state where the peak P1 is extracted. FIG. 11D is a diagram illustrating a state where the peak P2 is extracted. If there is no peak, the process proceeds to step S55 in FIG. If there is a peak, the process proceeds to step S54, and the peripheral pixel value of the extracted peak point Pn is cleared to 0 (see FIGS. 11E and 11F).

  FIG. 12A is a diagram illustrating processing subsequent to FIG. Proceeding from step S54 to step S55, the barycentric point Pg (Xg, Yg) is calculated for the original image (see FIG. 12D from FIG. 12B and FIG. 12C).

  The offset and angle (γ1, γ2) from the center of gravity Pg (Xg, Yg) of each peak point (Px1, Py1), (Px2, Py2) are calculated.

Next, pre-processing of a rendering image will be described with reference to FIGS.
FIG. 13B shows an example of a left-direction rendered image. The image size is about W300 × H300 pixels (the image size may slightly increase or decrease because it may protrude from the frame depending on the case), and the gradation is 256 monochrome gradations. is there. On the other hand, the result of binarization with a fixed threshold value (for example, 25 is set for 0 to 255 values) is shown in FIG. 13C. The inside of the eyeball is 1 and the background is 0. It is done. At this time, if the gradation is not uniform depending on the imaging state of the original image, a dust region may remain due to binarization as shown in the lower right. The dust of several pixels can be avoided by applying smoothing in advance using a 3 × 3 averaging matrix or the like before the binarization process as described in Patent Document 8 above. For a dust area with a relatively large area, a process of extracting a closed area from the binary image based on the connectivity of 8 neighboring pixels is performed, and the area is below a predetermined threshold (for example, set to 200 pixels) The process of deleting the closed region is performed. As a result, the inside of the eyeball is 2 and the background is 0. Finally, a contour pixel is set at the boundary between the inside of the eyeball and the background. Specifically, when four or more pixels are background pixels among eight neighboring pixels of the pixels set inside the eyeball, the value of the pixel is changed to 1 as a contour pixel. The result is shown in FIG. 13D, where 2 is given for the inside of the eyeball, 1 is given for the outline, and 0 is given for the background.

More specifically, FIG. 13 is a diagram for explaining details of the binarization processing and the contour extraction processing (S32) of the projection image. First, in step S32-1, a smoothing filter process is performed on the projected grayscale image in either the horizontal direction or the vertical direction (see FIG. 13B).
Here, the source image data has 256 gradations, and Image (x, y) = 0 to 255 (x = 0,..., Xs-1; y = 0,..., Ys-1).

In step S32-2, binarization processing of the projection image is performed (see FIG. 13B).
Here, the binary image data is
ImageM (x, y) = {0: background, 1: inside}.

Next, in step S32-3, a neighboring pixel connection process and a region extraction process are performed.
Here, the multi-value image data is
ImageM (x, y) = {0: Background, 2 or more: Internal} (x = 0, ..., Xs-1; y = 0, ..., Ys-1).

Next, in step S32-4 in FIG. 14A, the neighboring region is integrated (see FIG. 14B).
Here, the multi-value image data is
ImageM (x, y) = {0: Background, 2 or more: Internal} (x = 0, ..., Xs-1; y = 0, ..., Ys-1).

Next, in step S32-5, a minute area region is deleted (see FIG. 14C), the region 3 is deleted.
Here, the multi-value image data is
ImageM (x, y) = {0: Background, 2: Internal} (x = 0, ..., Xs-1; y = 0, ..., Ys-1).

Next, in step S32-6, contour extraction processing is performed (see FIG. 14D).
Here, the ternary image data is
Image (x, y) = {0: background, 1: contour, 2: inside} (x = 0, ..., Xs-1; y = 0, ..., Ys-1). In this way, the projected image of the eyeball can be distinguished into three values: background, contour, and internal.

FIG. 15A is a diagram illustrating an outline of an automatic setting function of an eye axis.
As described above, the eye axis serves as a reference for calculating parameters such as the left-right symmetry and vertical symmetry of the rear pole, and is required to be set with high accuracy. However, as long as the shape of the eyeball changes greatly depending on the pathological condition, it is difficult to automatically set the eye axis based on uniform eyeball shape characteristics. Intensity myopia, we propose a method to automatically set the eye axis based on the symmetry of the front pole, assuming that the symmetry of the front pole is maintained even if the back pole side is generally asymmetric. If this assumption does not hold, such as when the front pole is also asymmetric, manual interaction type correction is applied to deal with each case individually.

For the image of FIG. 14 (d), the value determined for one and becomes the minimum value X ₁ of the X direction of the pixel and a maximum value X _2, the Y-direction pixel position at which each takes the minimum and maximum Y ₁ and Y ₂ . In addition, an average value of X coordinates of pixels having a value of 1 or more is Xg, and an average value of Y coordinates is Yg. Based on these coordinate values, as shown in FIG. 15A (1), (X ₁ , Y ₁ ) is the front pole edge point Pf, (X ₂ , Y ₂ ) is the back pole edge point Pb, (Xg, Yg) is defined as the barycentric point Pg. At this time, Pf-Pg-Pb corresponds to the eye axis, but the eye axis temporarily set at this stage does not have an appropriate direction unless the anterior and posterior poles are symmetrical like the normal eye. .

The most reliable point among the three points Pf, Pg, and Pb is Pg. Therefore, Pf is corrected based on Pg. Specifically, as shown in FIG. 15A (2), the front pole edge point Pf is slightly moved up and down on the contour line of the front pole within a range where the X coordinate of Pf does not approach Pg. At each moved position, an intermediate measurement point Pc is set on the straight line Pg-Pf, a straight line is extended from Pc in the vertical direction perpendicular to the straight line Pg-Pf, and the intersections with the contour lines are defined as P _{+ v} and P _{−. v} and calculates the difference between the distance _{Pc-P + v} and _{Pc-P -v.} The position of Pf when this difference takes the minimum value is determined as the front pole edge point. At this time, Pc is not a single place, but is moved within a predetermined range on the straight line Pg-Pf (for example, Pc is moved by a distance of 2 to 12 pixels in the Pf direction from Pg), and each calculated distance Pc-P _{+ v} And the position of Pf where the sum of the differences between Pc-P- _v takes the minimum value is searched. Subsequently, Pg is corrected based on Pf determined in reverse. Specifically, as shown in FIG. 15A (3), the center of gravity Pg is moved slightly up and down with the X coordinate Xg fixed while the Y coordinate of Pg does not approach the upper and lower contour lines. Similarly to the above, an intermediate measurement point Pc is set on the straight line Pf-Pg at each moved position, a straight line is extended from Pc in the vertical direction perpendicular to the straight line Pf-Pg, and the intersection with the contour line is set. and P _{+ v} and _{P -v,} calculates the difference between the distance _{Pc-P + v} and _{Pc-P -v.} The position of Pg when this difference takes the minimum value is determined as the barycentric point. Similarly, Pc is not a single place, but is moved within a predetermined range on the straight line Pf−Pg (for example, Pc is moved by a distance of 2 to 12 pixels in the Pg direction from Pf), and each calculated distance Pc−P _{+ v} And the position of Pg where the sum of the differences between Pc-P- _v takes the minimum value.

  Finally, Pb is corrected so that Pf-Pg-Pb approaches a straight line with reference to the determined Pf and Pg. Specifically, as shown in FIG. 15A (4), the rear pole side edge point Pb is slightly moved up and down on the contour line of the rear pole within a range where the X coordinate of Pb does not approach Pg. At each moved position, the distance between the intersection point of the perpendicular line drawn from the barycentric point Pg to the straight line Pf-Pb and the straight line Pf-Pb and the barycentric point Pg is calculated. The position of Pb when this distance takes the minimum value is determined as the rear pole edge point.

With the above processing, Pf-Pb can be determined as the eye axis, but the eye position is displayed on the screen for judgment by the interpretation doctor. If inappropriate, the two points Pf and Pb are corrected to the interactive form on the screen. Then, the midpoint of the corrected Pf and Pb is reset to Pg. In addition, in order to perform parameter analysis, as shown in FIG. 7A (b), a measurement origin Po is set at a predetermined distance (for example, 87 pixels) from Pb on the determined eye axis. Intersections P _{+ w} and P _−w and P _{+ n} and P _−n on the rear polar contour of a straight line extending in two directions of angles (for example, ± θ _w and ± θ _n ) up and down around the edge point Pb Set.

More detailed processing contents will be described below.
FIG. 15B is a diagram illustrating details of the automatic setting process of the eye axis (step S33).
As shown in FIG. 15B (a1), in step S33-1, the edge points (X1, Y1), (X2, Y2) and the center of gravity ( Xg, Yg) is calculated. In the case of a projected image in the left-right direction (FIG. 15B (a2)), the maximum and minimum edge points in the X-axis direction are the front and rear pole side edge points, and in the case of the projected image in the vertical direction, the maximum and minimum edge points in the Y-axis direction are The edge point is on the pole side (step S33-2). FIG. 15B (b) to FIG. 15B (e) are diagrams showing a state in which processing related to the projected images in the horizontal and vertical directions is performed.

  FIG. 16A is a diagram showing the flow of processing following FIG. Here, it is assumed that the symmetry is maintained on the front pole side. In step S33-3, correction and determination processing of the front pole side edge point is performed (see FIG. 16B). Here, the front pole side edge point is slightly moved on the contour line, and the two distances (double arrows) extending from the center of gravity Pg to the contour line in the vertical direction from the vector in the front pole side edge point direction are closest. Thus, the front pole side edge point Pf is determined.

Next, in step S33-4, the center of gravity point is corrected and determined (see FIG. 16C).
Two distances (double arrows) in which the center of gravity Pg is moved up and down in the vertical direction inside the contour line and extended from the previously determined front pole side edge point Pf to the contour line in the vertical direction from the vector in the center of gravity point direction. The center of gravity Pg ′ is determined so that is closest.

  Next, as shown in FIG. 17A, in step S33-5, the rear pole edge point correction and the eye axis determination processing are performed (see FIG. 17B). Here, the rear pole side edge point Pb is determined so that the rear pole side edge point is slightly moved on the contour line and is closest to the extension line of the vector of the center of gravity Pg ′ from the predetermined front pole side edge point Pf. To do.

  FIG. 18A is a diagram showing a flow of eye axis correction processing. First, in step S35-1, either the front pole side edge point Pf or the rear pole side edge point Pb is selected on the screen (see FIG. 18B).

  Next, in step S35-2, the selected front pole side edge point Pf or rear pole side edge point Pb is moved to a desired position (Px) on the screen (see FIG. 18C).

  Next, automatic correction of the edge point Px set by movement onto the contour line LR and automatic correction of the center of gravity Pg ′ are performed (see FIG. 18D). The edge point is moved to the position on the nearest contour line, and the barycentric point is moved on the straight line connecting the updated edge point and the other non-updated edge point).

  In step S35-4, it is determined whether or not the direction is appropriate. If YES, the process ends. If NO, the process returns to step S35-1.

FIG. 19A is a diagram showing details of the measurement edge point creation process (S38) in the rear pole contour portion.
In step S38-1, 90 ° and −90 ° measurement edge points (Xp9, Yp9), (Xm9,...) On the upper and lower contour lines LR in the direction orthogonal to the eye axis C from the measurement center point (X0, Y0). Ym9) is added (see FIG. 19B).

In step S38-2, 90 °, 0 °, on the basis of -90 degree measuring edge points, theta _w and - [theta] _w of the measurement edge point (Xp4, YP4), to add (Xm4, Ym4) ( FIG. 19 (c)).

In step S38-3, θ _n and −θ _n degree measurement edge points (Xp2, Yp2) and (Xm2, Ym2) are added based on the θ _w , 0 degree, and −θ _w measurement edge points (FIG. 19 (d)).
In this way, measurement edge points are created at the rear pole.

  FIG. 20A is a flowchart showing the flow of the smoothing filter process (FIG. 11A: S52). FIG. 20B is a diagram illustrating an example of the smoothing filter.

  The filter operation is performed by multiplying the matrix with image data. A plurality of areas are selected and referred to near the pixel of interest to be processed (for example, 3 × 3, up, down, left, and right). This value (1/9) is multiplied by the pixel value, and an average value is obtained by the addition process (an average value of the sum total of 9 pixels of 3 × 3).

  As shown in FIG. 20B, a 3 × 3 smoothing filter matrix is shown as an example of M (i, j). In step S52-1, a 3 × 3 smoothing matrix shown in FIG. 20B is defined. Next, in step S52-2, initialization is defined as x = y = 1, i = j = -1, and a temporary variable (temporary variable of pixel value) v = 0. v is x, and y is an x-direction address and y-direction address (total Xs, Ys pixels) for scanning pixels. i and j are displacements to be applied to calculate the matrix, the input projection image is Image (x, y), x = (0, 1,..., Xs−1), y = (0, 1, ..., Ys-1). i and j take integer values ranging from −1 to +1.

  In step S52-3, v ← v + Image (x + i, y + j) · M (i, j) is set. In step S52-4, i is incremented by 1. In step S52-5, i and +1 are compared. If i> 1, the processing by one 3 × 3 filter is finished, and step S52- Proceeding to step 6, j is incremented by 1, and the same processing is performed to finish the processing of multiplying the 3 × 3 matrix. Then, in Steps S52-9 and S52-11, it is determined whether or not calculation has been performed for all pixels, and the process ends (End). In this manner, by moving the 3 × 3 smoothing matrix on the input projection image, the pixel values of the projection image can be smoothed by the 3 × 3 smoothing filter.

  FIG. 21 is a flowchart showing the flow of the shading peak point extraction process (see FIG. 11A: S53) for the projected image. The point that loop processing is performed with x and y is the same as in FIG. However, processing up to Xs and Ys is performed. Pixel Image (x, y) has a value from 0 to 255. This is a process of finding a local peak having a local peak in the projection image. The pixel has a value from 0 to 255. In step S53-1, as an example of the initial value, Vm = 128 is set as a value not smaller than 128, and the initial value (-1, x) is set as the peak coordinates (Px, Py). -1) is set. Further, a temporary copy Image2 (x, y) obtained by copying the grayscale image Image (x, y) is provided as a backup, and x = y = 0 is set as an initial value. In step S53-2, Image2 (x, y) is compared with the light and dark peak Vm. If Image2 (x, y)> Vm, the process proceeds to step S53-3, Vm = Image2 (x, y) is set, and Px = x and Py = y are saved. Hereinafter, when processing is performed for all pixels, a peak point is obtained. That is, if Image2 (x, y) <= Vm, step S53-3 is passed, and x is incremented by 1 in step S53-4. The process is repeated until the end of the image in the horizontal direction (x = Xs) (step S53-5), and then y is also processed (steps S53-6, 7). In step S53-8, it is determined whether or not Xp <0 and Yp <0. If Yes, it is determined that there is no dark peak (no peak point), and the initial value is reset, for example. If No, there is a peak (Xp, Yp), and this peak (Xp, Yp) is output.

Next, the process proceeds to FIG. FIG. 22 is a flowchart showing the flow of processing for clearing the pixels around the peak point extracted by the processing of FIG. 21 (step S54 of FIG. 11). From the peak output in FIG. 21, in step S54-1, the peak coordinates (Xp, Yp) are reset as x = y = 0. In step S54-2, the distance r to the peak point (Xp, Yp) is set to r = {(x−Px) ² + (y−Py) ² } ^1/2.
Calculated by In step S54-3, the distance r is compared with a threshold value (for example, 50 pixels). If r is equal to or smaller than Sr, the process proceeds to step S54-4, where Image2 (x, y) = 0, and when r is larger than Sr, x is incremented by 1 in step S54-5, and in step S54-6 , X and Xs are compared. If x <Xs, the process returns to step S54-2. If x is equal to or greater than Xs, the process proceeds to step S54-7, the same process is performed for y, and the clear process is terminated. Thereby, it is possible to clear the pixel values in the region in the range from the peak point (Xp, Yp) to the radius Sr.

  FIG. 23 is a flowchart showing the calculation process of the barycentric point for the original image (FIG. 12A: step S55). Here, the barycentric point is obtained by taking the average value of the pixels.

First, in step S55-1, the center of gravity coordinates (Xg, Yg) are initialized as Xg = Yg = 0, weight sum V = 0, x = y = 0. Next, in step S55-2,
Xg ← Xg + Image (x, y) ・ x, Yg ← Yg + Image (x, y) ・ y, V ← V + Image (x, y)
As a result, the average value of the pixels is calculated by the processing from step S55-3 to step S55-6, and in step S55-7,
Xg ← Xg / V, Yg ← Yg / V
Then, the barycentric point for the original image can be obtained (see FIG. 12D).

  FIG. 24 is a flowchart showing details of pixel binarization processing (step S32). Pixel values 0 to 255 are converted to either 0 or 1. That is, if it is larger than Sl, it is set to “1”, otherwise it is set to “0”.

  First, in step S32-2-1, binary image data ImageM (x, y) = 0 (initial value is all pixels 0) and x = y = 1 are set. In step S32-2-2, Image (x, y) is compared with a threshold value Sl (L). The threshold value Sl is 25, for example. If Image (x, y) is larger than the threshold value Sl, the process proceeds to step S32-2-3, and binary image data ImageM (x, y) = 1 is set, and Image (x, y) is the threshold value. When it is equal to or less than Sl, the process proceeds to step S32-2-4, x is incremented by 1, and the above process is repeated until x becomes Xs-1 as shown in step S32-2-5. . When x becomes larger than Xs-1, the process proceeds to step S32-2-6, y is incremented by 1, and the process is repeated until y becomes larger than Ys-1. When y becomes Ys-1 or more, the process ends. Thus, it is possible to perform binarization processing of an image that is set to “1” if it is larger than S1 (L), and to “0” otherwise (see FIG. 13C).

  FIG. 25A is a flowchart showing the flow of neighboring pixel connection processing and region extraction processing (see step S32-3, FIG. 13D). FIG. 25B is a diagram illustrating four neighboring pixels among the eight neighboring pixels. When binarization is performed with 0 and 1, for example, the value of the target pixel is 1 and whether or not a pixel having a value (region Id) of 2 or more exists in any of the four neighboring pixels of the target pixel . When there is a neighboring pixel having a value of 2 or more, the value of the target pixel is set to the Id of the neighboring pixel, and the area value (number of pixels in the region) in the region table corresponding to the Id of the neighboring pixel Is incremented by one. When there is no neighboring pixel having a value of 2 or more, a new area Id (initial value is 2) is assigned to the target pixel, a new record is registered in the area table corresponding to the Id, and the initial area value is set. Set to 1. By repeating such processing, two or more areas Id are allocated to all pixels having 1 in the image binarized by 0 and 1, and are divided into N areas, and are also stored in the area table. The area calculation result of each area is stored.

  First, in step S32-3-1, Id: region Id (initial value: 1), region table Tb (id): initial value 0, n: number of regions (initial value n = 0), x = y = 1 Set. In step S32-3-2, whether or not the target pixel (x, y) is inside the eyeball, that is, ImageM (x, y) is compared with 0. If it is inside the eyeball, that is, if it is larger than 0, in step S32-3-3, four 8-neighbor pixels (x-1, y-1), (x, y-1), (x + 1) , y-1), (x-1, y). That is, for the four pixels whose addresses are smaller than their own addresses, the background (0) and the inside of the eyeball (2 or more, there is no value 1 for processed pixels with a small address) are designated as internal. Check if it exists. That is, in step S32-3-4, it is determined whether there is a pixel having a value Idp of 2 or more.

In the case of Yes, the Idp of the neighboring pixel is diverted as the area Id of its own pixel. In the case of No, the process proceeds to Step S32-3-6 to assign a new area Id, and Id ← Id + 1, Idp = Id, n ← n + 1, and in the case of Yes, the process proceeds to step S32-3-5, and 1 is added to the count value (initial value 0) of the pixel to the table value (No. 2), and at the same time The image data is updated as the area number Idp (= 2).
Tb (Idp) ← Tb (Idp) +1, ImageM (x, y) = Idp

  Next, the process proceeds to step S32-3-7, and x is incremented by 1. In step S32-3-8, x and Xs-1 are compared, and the process is repeated until the contour, and the same process is performed for y. I do. Thus, the neighboring pixel connection processing and region extraction processing are terminated (End: see FIG. 13D). Thereby, area division and area calculation of each area are completed. However, since only four neighboring pixels are referenced when assigning the region Id, different regions Id are assigned to neighboring regions to which the same Id should be originally assigned, and the region tends to be excessively divided. Therefore, correction processing for integrating neighboring regions is performed.

  FIG. 26A is a flowchart showing the flow of the neighborhood region integration process (see FIG. 14, step S32-4, FIG. 14B). FIG. 26B is a diagram illustrating an example of 8 neighboring pixels. This is a process of organizing the extracted regions by looking at all eight neighboring pixels in a state where the region table is completed, and the region numbers of all the pixels are reassigned and organized.

As a result of the above processing, the multi-value image data is
ImageM (x, y) = {0: Background, 2 or more: Internal} (x = 0, ..., Xs; y = 0, ..., Ys)
The pixels having values of 2 or more are integrated so as to be as small as possible (within a range of 2 or more).

  First, in step S32-4-1, x = y = 1, and in step S32-4-2, whether or not the target pixel (x, y) is inside the region, that is, ImageM (x, y) is determined. Determine whether it is positive or zero. If it is positive, the process proceeds to step S32-4-3, and 8 neighboring pixels (x-1, y-1), (x, y-1), (x + 1, y-1), (x-1 , y), (x + 1, y), (x-1, y + 1), (x, y + 1), (x + 1, y + 1). That is, in step S32-4-4, it is determined whether there are pixels in different regions in the vicinity where ImageM (x + i, y + j) <ImageM (x, y). That is, a pixel in the vicinity that has a value of 2 or more and is numbered lower than its own ID number is examined. In the case of Yes, the process proceeds to step S32-4-5, and the area Id update process is performed. In the case of No, in the same manner as in the case where ImageM (x, y) is 0 in step S32-4-2, Proceed to 4-6, x is incremented by 1, and the process is continued. When x becomes Xs−1, the same processing is performed with respect to y, and when y = Ys−1, the processing ends (End).

  Next, as shown in FIG. 27, when a pixel in the vicinity that is assigned a younger number than its own ID number is examined and a young number is assigned, its own pixel is re-entered into the younger number (younger) Group). When the own area id is replaced with the Idp of the peripheral pixel (x + i, y + j), all the pixels having the same area id as the self are changed to Idp over the entire image data. Then, since there is no pixel having the value of the region id, the values of all the pixels having a value larger than the value of the region id are decreased by 1. As a result, the type N of the area id to be used is decreased by 1, and by repeating this process, the areas are sequentially integrated and N becomes smaller.

In region step S32-4-11,
Set id = ImageM (x, y), Idp = ImageM (x + i, y + j), u = v = 1, and compare ImageM (u, v) with id in step S32-4-12 To do. If ImageM (u, v) is 0, the process proceeds to step S32-4-15. If ImageM (u, v) is id, the process proceeds to step S32-4-13, where ImageM (u, v) = Idp, and the process proceeds to step S32-4-15. If ImageM (u, v)> id, the process proceeds to step S32-4-14, and ImageM (u, v) ← ImageM (u, v) -1 is set, and the process proceeds to step S32-4-15. In step S32-4-15, u is incremented by 1. In step S32-4-15, u and Xs-1 are compared. If u ≧ Xs-1, the process proceeds to step S32-4-17. If u <Xs-1, the process returns to step S32-4-12. In step S32-4-17, u = 1, v ← v + 1, and v and Ys-1 are compared in step S32-4-18. If v ≧ Ys−1, the process proceeds to step S32-4-19. If v ≧ Ys−1, the process returns to step S32-4-12. In step S32-4-19, regarding the area table Tb, the area id is integrated into Idp, and the id is cleared to 0. That is,
Tb (Idp) ← Tb (Idp) + Tb (id), Tb (id) = 0, id ← id + 1
In step S32-4-20, id and n + 2 are compared. If id ≧ n + 2, the process proceeds to step S32-4-22. If id <n + 2, step S32- Proceed to 4-21. In step S32-4-21, the id numbers of all areas are decreased by 1. That is,
Tb (id-1) ← Tb (id), id ← id + 1
Return to Step S32-4-20. When id ≧ n + 2 is finally reached, n is decremented by 1 in step S32-4-22, and the process is terminated. As a result, it is possible to perform integration processing of neighboring regions as shown in FIG.

  FIG. 28 is a flowchart showing a deletion process of a small area region (step S32-5) for deleting a small area such as peripheral dust and reducing an error in determining an eye axis (FIG. 14). (See (c)). First, id = 2 is set in step S32-5-1, and Tb (id) and St are compared in step S32-5-2. Here, St is a lower limit threshold value of the area, for example, 200. If Tb (id)> = St, the process proceeds to step S32-5-10. If Tb (id) <St, the process proceeds to step S32-5-3, and all pixels having the value id in the image data are set as background pixels. In step S32-5-3, the variables u and v are set to 1. In step S32-5-4, ImageM (u, v) is compared with the area number id. If ImageM (u, v) = id, in step S32-5-5, ImageM (u, v) = 0 and the pixels in the area are set as the background area, and the process proceeds to step S32-5-6. If ImageM (u, v) is not equal to id, nothing is done and the process proceeds to step S32-5-6.

  In step 32-5-6, u is incremented by 1. In step 32-5-7, u is compared with Xs-1. If u <Xs−1, the process returns to step S32-5-4. If u> = Xs−1, in step S32-5-8, u = 1 and v ← v + 1.

  Next, in step S32-5-9, v and Ys-1 are compared. If v≥Ys-1, the process returns to step S32-5-4, and if v≥Ys-1, step S32-5-5 is performed. Proceeding to 10, Id is incremented by 1, and in step S32-5-11, id is compared with n + 2. If id <n + 2, the process returns to step S32-5-4, and if id ≧ n + 2, the process proceeds to step S32-5-12.

  As described above, after a region having a small area is deleted in the previous processing, if a value of region number 3 or more still remains, the value of 3 or more is unconditionally changed to 2 (trinarization is performed, FIG. 14). (See (c)). That is, in step S32-5-12, u = v = 1, and in step S32-5-13, the size of ImageM (u, v) is compared with 2, and if ImageM (u, v)> 2, In step S32-5-14, ImageM (u, v) = 2 is set, and in the case of ImageM (u, v) <= 2, u is incremented by 1 in step S32-5-15. As a result, the three or more areas in FIG. In step S32-5-16, u is compared with Xs-1, and if u <Xs-1, the process returns to step S32-5-13, and if u> = Xs-1, step S32-5 is returned. Proceed to 17, and set u = 1 and v ← v + 1. Next, in step S32-5-18, v and Ys-1 are compared. If v <Ys-1, the process returns to step S32-5-13, and if v> = Ys-1, processing is performed. finish. As a result, an area smaller than a certain threshold value St is deleted (the process of replacing the area No. 3 in FIG. 14B with the background (0)), that is, an area having a very small area can be deleted.

  FIG. 29 is a flowchart showing the flow of the contour extraction process (see FIG. 14D). A process of counting the number of pixels that are 2 (internal) in a total of 9 pixels in the vicinity of 8 including myself, and changing it to “1” as an outline if it is a value from 3 to 7 It is. In step S32-6-1, x = y = 1 is set. Next, c = 0 and i = j = −1 are set in the pixel counter. In step S32-6-3, whether or not eight neighboring pixels (x + i, y + j) of the target pixel (x, y) are inside the region, that is, ImageM (x + i, y + j) and If ImageM (x + i, y + j)> 0 is compared with 0, c is incremented by 1 in step S32-6-4, and in the case of = 0 in step S32-6-3 In step S32-6-5, i is incremented by one. Next, i is compared with 1, and if i ≦ 1, the process returns to step S32-6-3. If i> 1 is compared with i, the process proceeds to step S32-6-6 and i is incremented by 1. In step S32-6-6, i and +1 are compared. If i ≦ 1, the process returns to step S32-6-3. If i> 1, i = −1 in step S32-6-7. , j ← j + 1.

  Next, in step S32-6-8, j and +1 are compared. If j ≦ 1, the process returns to step S32-6-3. If j> 1, the process returns to step 32-6-9. It is determined whether the target pixel has a contour based on the value of the number c of neighboring pixels in the counted area. For example, it is determined whether or not 3 ≦ c ≦ 7 (the contour is determined based on whether the number of 8 neighboring pixels falls within the range of 3 to 7. In the case of yes, in Step S32-6-10, ImageM ( x, y) = 1, and when it is no, in step S32-6-11, x is incremented by 1. Then, in step S32-6-12, x and Xs-1 are compared. If x <Xs-1, the process returns to step S32-6-2, and if x ≧ Xs-1, the process proceeds to step S32-6-13, where x = 1, y ← y + 1, and step S32- 6-14, y and Ys-1 are compared, and if y≥Ys-1, the process returns to step S32-6-2, and if y≥Ys-1, the process is terminated to extract the contour. can do.

  FIG. 30 is a flowchart showing a flow of processing for calculating the front pole side and rear pole side edge points and the barycentric point (S33-1: FIG. 15B). Although it is the same as the process for obtaining the center of gravity described above, it takes only three values 0, 1, and 2 instead of a gradation image from 0 to 255.

  First, in step S33-1-1, barycentric coordinates (Xg, Yg), Xg = Yg = 0, weight sum V = 0, x = y = 0, front pole side edge point (X1, Y1), X1 = Y1 = -1, Set the rear pole edge point (X2, Y2), X2 = Y2 = -1.

  Next, in step S33-1-2, whether or not the target pixel (x, y) is inside the region or in the outline, that is, ImageM (x, y) is compared with 0. If ImageM (x, y) = 0, the process proceeds to step S33-1-13. If ImageM (x, y)> 0, the process proceeds to step S33-1-3, where Xg ← Xg + x, Yg ← Yg + y, V ← V + 1. Next, in step S33-1-4, it is determined whether or not the projected image is in the left-right direction. In the case of Yes, it progresses to step S33-1-5, and in No, it progresses to step S33-1-9.

  In step S33-1-5, it is determined whether or not X1 <0 or x <X1. If Yes, X1 = x, Y1 = y (minimum value candidate) is set in step S33-1-6. In the case of No, the process proceeds to step S33-1-7 to determine whether X> X2, and in the case of Yes, X2 = x, Y2 = y. Along with No, X2 = x, Y2 = y (maximum value candidate). In this way, the process is to find the minimum and maximum in the x direction.

In the case of No in step S33-1-4, similarly in the vertical direction, the process proceeds to step S33-1-9 to determine whether Y1 <0 or y <Y1, and in the case of Yes, Proceeding to step S33-1-19, X1 = x, Y1 = y is set, and in addition to the case of No, in step S33-1-11, it is determined whether y> Y2. In the case of Yes, X2 = x and Y2 = y are set in step S33-1-12, and with the case of No, 1 is incremented to x in step S33-1-13. In step S33-1-14, x and Xs are compared. If x <Xs, the process returns to step S33-1-2, and if x ≧ Xs, x = 0, y in step S33-1-15. ← y + 1, and in step S33-1-16, y is compared with Ys. If y ≧ Ys, the process returns to step S33-1-2. If y ≧ Ys, in step S33-1-17, Xg ← Xg / V, Yg ← Yg / V is set, and step S33-1-18. In step S33-1-19, X1 << X2, Y1 << is determined. In step S33-1-19, the process is terminated. The process ends as Y2. Thereby, X1, X2, Y1, and Y2 can be obtained.
Corrections are made in order to the temporarily determined points.

  FIG. 31 is a flowchart showing a flow of correction and determination processing of the front pole side edge point (see S33-3 and FIG. 16). In step S33-3-1, in order to correct X1 and Y1, a distance difference d is calculated from (X1, Y1) and (Xg, Yg) by a method described later. First, in step S33-3-2, a front pole side edge point (x1, y1), an initial value x1 = X1, y1 = Y1, and a minimum distance difference dm = d, x = y = 0 are set.

In step S33-3-3, whether or not the target pixel (x, y) is on the contour, that is, ImageM (x, y) is compared with 1, and if not, it is determined that it is not a point on the contour. Proceed to step S33-3-8, and if it is the same, determine that the point is on the contour, and proceed to step S33-3-4,
Left direction: x <{X1 * (Sw-1) + Xg} / Sw
Right direction: x> {X1 * (Sw-1) + Xg} / Sw
Up direction: y <{Y1 * (Sw-1) + Yg} / Sw
Down: y> {Y1 * (Sw-1) + Yg} / Sw
It is determined whether or not the above is satisfied. Here, Sw is a scanning weight, for example, 16. Here, a range of points in the candidate is determined.

  Here, if it is No, it will progress to step S33-3-8. If Yes, in step S33-3-5, the distance difference d is calculated from (x, y) and (Xg, Yg) (see FIG. 16B). Here, in step S33-3-6, it is determined whether d <dm. In the case of Yes, it progresses to step S33-3-7, and it updates with x1 = x, y1 = y, dm = d (front pole side edge point candidate), and in the case of No, step S33-3-8 Proceed to By repeating such processing, the edge point on the front pole side is determined.

  In step S33-3-8, x is incremented by one. Next, in step S33-3-9, x and Xs are compared. If x <Xs, the process returns to step S33-3-3. If x ≧ Xs, the process proceeds to step S33-3-10, y and Ys are compared, and if y ≧ Ys, the process returns to step S33-3-3. If y ≧ Ys, the process proceeds to step S33-3-12, and X1 ← x1, Y1 ← y1, and the correction processing for the front pole side edge point is completed.

FIG. 32 is a flowchart showing a flow of a distance difference calculation process between the front pole side edge point and the gravity center point, which is the process of steps S33-3-1 and S33-3-5 of FIG. 31 (FIG. 16, step S33-4). In FIG. 32, since this process is used not only for correcting the front pole side edge point but also for correcting the center of gravity point described later, the front edge point with respect to the front pole side edge point and the rear edge with respect to the center of gravity point. It has a generalized notation called point. Both (X1, Y1) in step S33-3-1 and (x, y) in S33-3-5 in FIG. 31 correspond to the front edge point (Xs, Ys) in FIG. Both (Xg, Yg) in S33-3-1 and S33-3-5 correspond to the rear edge point (Xe, Ye) in FIG. First, in step S33-4′-1, the distance difference is calculated.
The front edge point (Xs, Ys), the rear edge point (Xe, Ye),
Unit vector in the direction from the front edge point to the rear edge point
dx = (Xe-Xs) / L and dy = (Ye-Ys) / L are obtained.

However, L = {(Xe-Xs) ² + (Ye-Ys) ² } ^1/2
Define distance r (initial value: 2, upper limit: 12) from front pole side edge point P _f ,
The initial value = 0 of the output distance difference is calculated 11 times from r = 2 to = 12.

In step S33-4′-2,
Specify the measurement intermediate point xc = dx × r + Xs, yc = dy × r + Ys (the measurement intermediate point (xc, yc) specified on the axis from the front edge point to the rear edge point is obtained .), Extending a straight line upward from the measurement intermediate point perpendicular to the axis from the front edge point to the rear edge point, and calculating the intersection (xc1, yc1) on the upper contour of the eyeball region as follows: .

In step S33-4'-3, variables are initialized by substituting xc1 = yc1 = -1, t = 0.5 (0.5 pixel unit), and in step S33-4'-4,
x = -dy × t + xc, y = dx × t + yc, that is, the unit vector in the direction orthogonal to the axis from the front edge point to the rear edge point is changed by the value t Move up. While increasing t, the line is extended upward to search for a place where the outline collides, that is, as shown in step S33-4′-5,
It is found by determining whether x <Xs, y <Ys, and ImageM (x, y) = 1.

In the case of Yes, in step S33-4′-6, the point where the line extending in the upward direction and the outline intersect (FIG. 16B).
xc1 = x, yc1 = y,
In the case of No, 0.5 is added to t in step S33-4′-7, and the process returns to step S33-4′-4.

Next, a straight line is extended in the downward direction perpendicular to the axis from the front edge point to the rear edge point from the measurement intermediate point, and the intersection point (xc2, yc2) on the lower contour of the eyeball region is similarly calculated. The process proceeds from step S33-4′-6 to step S33-4′-8, where xc2 = yc2 = −1 and t = −0.5. Next, in step S33-4′-9,
x = −dy × t + xc, y = dx × t + yc.

Next, in step S33-4′-10,
It is determined whether x <Xs, y <Ys and whether ImageM (x, y) = 1. If No, 0.5 is subtracted from t in step S33-4′-12, and step S33. Return to -4'-9. In the case of Yes, in step S33-4′-11, the point where the extended line and the outline intersect (FIG. 16B).
xc2 = x, yc2 = y
In step S33-4′-13, the distance difference between the distance to the intersection on the contour above the measurement intermediate point and the distance to the intersection on the contour below the measurement intermediate point is expressed as follows: And is added to the sum d of distance differences.
d ← d + | (xc1-xc) ² + (yc1-yc) ^2- (xc2-xc) ^2- (yc2-yc) ² |
(Here, the obtained d is the sum of the distances between (xc1, yc1) and (xc2, yc2) with reference to the point (xc, yc)), and step S33-4 ′. −14, r is moved in the center direction (center of gravity), and in step S33-4′-15
r <Upper limit value (12)
By determining whether or not, the distance difference in the vertical direction is obtained over a plurality of locations while moving the measurement intermediate point (xc, yc), and the sum is d, and d is output.

  If yes, the process returns to step S33-4'-2, and if no, the process ends (End). By this processing, the distance difference d between the front pole side edge point and the center of gravity can be calculated.

  FIG. 33 is a flowchart showing the flow of the center-of-gravity point correction and determination process (S33-4) (see FIG. 16C). The basic idea is the same as that in FIG. 31. In FIG. 31, the center of gravity is provisionally determined, and the front pole side edge point is obtained from the provisional center of gravity. The time when the distance difference d becomes the smallest is the front pole side edge point. The difference is that in FIG. 31 it passes on the contour line, whereas the center of gravity varies internally. In step S33-4-1, the distance difference d is calculated from (Xg, Yg) and (X1, Y1). Each of (Xg, Yg) corresponds to the front edge point (Xs, Ys) in FIG. 32, and (X1, Y1) corresponds to the rear edge point (Xe, Ye) in FIG. The difference d is calculated. Here, in step S33-4-2, the barycentric point (xg, yg), the initial value xg = Xg, yg = Yg, and the minimum distance difference dm = d are set. Next, in step S33-4-3, it is determined whether or not the direction is the left-right direction. If No (up and down direction), the process proceeds to step S33-4-11. If Yes (left and right direction), the process proceeds to step S33-4-4.

In step S33-4-4 (left-right direction), the center of gravity is moved in the y-axis direction (the x-axis direction is not moved (xg is not changed with the initial value)).
Ymin = {Yg * (100-Sw) + Y1 * Sw} / 100
Ymax = {Yg * (100-Sw) + Y2 * Sw} / 100, y = Ymin
And As a result, a range to be moved is determined so that an unnecessary range is not moved. Here, Sw: the scanning weight is 30, for example.

  Next, in step S33-4-5, the distance difference d is calculated from (Xg, y) and (X1, Y1) in the Y direction. (Xg, y) corresponds to the front edge point (Xs, Ys) in FIG. 32, (X1, Y1) corresponds to the rear edge point (Xe, Ye) in FIG. The difference d is calculated.

First, in step S33-4-6, it is determined whether d <dm. In the case of Yes, in step S33-4-7,
yg = y, dm = d (center of gravity candidate)
And in the case of No, in step S33-4-8
y ← y + 1
In step S33-4-9, y> Ymax is determined.

If No, the process returns to Step S33-4-6. If Yes, the process returns to Step S33-4-10.
Yg = yg.

On the other hand, in step S33-4-11, Sw is a scanning weight (example 30). As a result, the range of movement in the vertical direction is determined so that the unnecessary range is not moved.
Xmin = {Xg * (100-Sw) + X1 * Sw} / 100
Xmax = {Xg * (100-Sw) + X2 * Sw} / 100, x = Xmin

In step S33-4-12,
For the X direction, the distance difference d is calculated from (x, Yg) and (X1, Y1). (x, Yg) corresponds to the front edge point (Xs, Ys) in FIG. 32, (X1, Y1) corresponds to the rear edge point (Xe, Ye) in FIG. The difference d is calculated.

Next, in step S33-4-13, d <dm is determined. In the case of No, the process returns to Step S33-4-13. In the case of Yes, the process proceeds to Step S33-4-17, and Xg is obtained by setting Xg = xg.
Thereby, as shown in FIG.16 (c), the correction | amendment and determination process of a gravity center point are complete | finished.

Next, FIG. 34 is a flowchart showing the flow of correction of the posterior pole side edge point and the determination process of the eye axis (automatic setting process S33 of the eye axis in FIG. 17). Since Pf and Pg ′ are determined, Pb is then obtained so as to be on an extension line of Pf and Pg ′. First, in step S35-5-1, the linearity d is determined from (X1, Y1)-(Xg, Yg)-(X2, Y2). In step S35-5-2, a rear pole edge point (x2, y2), an initial value x2 = X2, y2 = Y2, a minimum judgment value dm = d, and x = y = 0 are set. In step S33-5-3, ImageM (x, y) is compared with 1. If ImageM (x, y) is not 1, the process proceeds to step S33-5-8. If ImageM (x, y) is 1, the process proceeds to step S33-5-4. Here, it is confirmed whether or not the values of x and y are included in a predetermined range indicated by the following conditions.
Left direction: x> {X2 * (Sw-1) + Xg} / Sw
Right direction: x <{X2 * (Sw-1) + Xg} / Sw
Upward: y> {Y2 * (Sw-1) + Yg} / Sw
Downward: y <{Y2 * (Sw-1) + Yg} / Sw
Meet?
(Sw: scanning weight (example is 4)).

  In No, it progresses to step S33-5-8. In the case of Yes, in step S33-5-5, a process of determining the linearity d from (X1, Y1)-(Xg, Yg)-(x, y) is performed.

First, in order to determine a constraint, d <dm is determined in step S33-5-6. In No, it progresses to step S33-5-8. In the case of Yes,
x2 = x, y2 = y, dm = d (rear pole edge candidate)
It becomes. Next, in step S33-5-8,
x ← x + 1
In step S33-5-9, x is compared with Xs. If x <Xs, the process returns to step S33-5-3. If x ≧ Xs,
As x = 0, y ← y + 1, y and Ys are compared in step S33-5-11. If y ≧ Ys, the process returns to step S33-5-3, and if y ≧ Ys, in step S33-5-12,
X2 ← x2, Y2 ← y2
Then, the process ends (End). Thereby, X2 and Y2 can be corrected (see FIG. 17B).

  FIG. 35 is a flowchart showing the flow of the process of determining the linearity of the three vertices, which is a specific example of the processing contents of S33-5-1 and S33-5-5 of FIG. 34 (FIG. 17B). )reference). In FIG. 34, (X1, Y1) in S33-5-1 is (xs, ys) in FIG. 35, (Xg, Yg) is in (xm, ym) in FIG. 35, and (X2, Y2) is in FIG. (xe, ye), (X1, Y1) in S33-5-5 is (xs, ys) in FIG. 35, (Xg, Yg) is (xm, ym) in FIG. y) corresponds to (xe, ye) in FIG. First, in step S33-5-1-1, given three vertices are set to (xs, ys), (xm, ym), (xe, ye). Also, dx = xe-xs, dy = ye-ys, and output determination value d = 0.

  In step S33-5-1-2, | dx |> | dy | is determined. If No (up and down direction), the process proceeds to step S33-5-1-4. If Yes (left and right direction), the process proceeds to step S33-5-3.

In step S33-5-1-3, d is obtained by d = | dy × (xm-1) / dx + ys-ym |. In step S33-5-1-4, d = | dx × (ym -1) d is obtained from / dy + xs-xm |, and the process is terminated (End: see FIG. 17B)).
The process up to this point is the eye axis setting process.

Thereafter, manual adjustment is performed.
FIG. 36 is a flowchart showing an automatic correction process (S35-3) for manually setting edge points on the contour line. This is a process of finding a point on the contour line in FIG.

First, in step S35-3-1,
The front pole edge point to be corrected is (x1, y1), the corrected front pole edge point is (X1, Y1), the center of gravity is (Xg, Yg), and the rear pole edge point (X2, Y2) Is fixed. Although specific explanation is omitted, conversely, the front pole side edge point (X1, Y1) is fixed, the rear pole side edge point to be corrected is (x2, y2), and the rear pole side edge after correction A method in which the point is (X2, Y2) and the center of gravity is (Xg, Yg) can be similarly realized. The unit vector in the axial direction is defined as follows using the front pole side edge point before correction.
L = {(x1-X2) ² + (y1-Y2) ² } ^1/2
dx = (x1-X2) / L, dy = (y1-Y2) / L
t = 0.5
The corrected front pole side edge point is (x, y), and the initial values are x = x1 and y = y1.

In step S35-3-2, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether or not ImageM (x, y) = 1. In the case of Yes, it progresses to step S35-3-10, and in the case of No, it progresses to step S35-3-2.
x = dx × t + x1, y = dy × t + y1
And change. Next, in step S35-3-4, whether or not the target pixel (x, y) is on the contour, that is,
ImageM (x, y) = 1
In the case of Yes, the process proceeds to step S35-3-10. In the case of No, the process proceeds to step S35-3-.
x = -dx × t + x1, y = -dy × t + y1
And change. Next, in step S35-3-7, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether or not ImageM (x, y) = 1. In the case of Yes, it progresses to step S35-3-10, and in the case of No, it progresses to step S35-3-8, and is set to t ← t + 0.5. Next, in step S35-3-9,
Determine if t <L / 2,
If No, the process returns to Step S35-3-3, and if Yes, Error (cannot be corrected) is set. In Step S35-3-10,
X1 = x, Y1 = y
Xg = (X1 + X2) / 2
Yg = (Y1 + Y2) / 2
Then, the process ends (End). Thereby, automatic correction processing is performed on the contour line of the manually set edge point.

FIG. 37 is a flowchart showing the flow of processing for adding 90-degree and -90-degree edge points in a direction orthogonal to the already determined eye axis (see step S38-1, FIG. 19A). First, in step S38-1-1, a unit vector on the eye axis in the rear pole direction with the front pole side edge point as a base point is obtained. 36. Xg and Yg up to FIG. 36 are the geometric center of gravity of the eyeball introduced to automatically set the eye axis, and Xo and Yo in FIGS. 37 and 38 are the center of gravity for measurement. Since the values are a certain distance away from the set rear pole on the eye axis, the values are different between the two (however, the center of gravity for measurement is set so that they are almost the same for normal eyes). The unit vector on the eye axis is
dx = (X2-X1) / L, dy = (Y2-Y1) / L
However, L = {(X2-X1) ² + (Y2-Y1) ² } ^1/2
And

In step S38-1-2, t = 0.5. In step S38-1-3, as a point extended from the measurement center of gravity in the direction of the unit vector,
x = -dy × t + X0, y = dx × t + Y0
Define

Next, in step S38-1-4, scanning is performed in the upward direction, and a point hitting the contour line is obtained. Whether the target pixel (x, y) is on the contour, ie,
x <Xs, y <Ys, ImageM (x, y) = 1
In the case of No, in step S38-1-5,
t¬t + 0.5
And return to Step S38-1-3.

In the case of Yes, it progresses to step S38-1-6,
Xp9 = x, Yp9 = y
And
In step S38-1-7, t = 0.5.

Next, in step S38-1-8, the scan is performed in the downward direction, and a point that hits the contour line is obtained.
x = -dy × t + X0, y = dx × t + Y0,
In step S38-1-9, whether or not the target pixel (x, y) is on the contour, that is,
It is determined whether x <Xs, y <Ys, and ImageM (x, y) = 1. In No, it progresses to step S38-1-10, and after setting it as t¬t-0.5, it returns to step S38-1-8.

In the case of Yes, in step S38-1-11,
Xm9 = x, Ym9 = y
And
Thereby, an edge point for measurement can be added.

Figure 38 is a flowchart showing the flow of a process of adding theta _w degree and - [theta] _w of the measuring edge point (± first angle) (Step S38-2, FIG 19 (a), (c) reference). Here, an example in which _θw is 45 degrees is shown, but the first angle is narrower than 90 degrees and wider than a second angle described later (40 degrees to 90 degrees).

First, in step S38-2-1 (measuring edge point plus 45 degrees), calculates the unit vector of theta _w degree direction. In this embodiment, dx and dy are unit vectors in the 45 degree direction.
dx = {(X2-X0) / L + (Xp9-X0) / Lp} / 2
dy = {(Y2-Y0) / L + (Yp9-Y0) / Lp} / 2
However, L = {(X2-X0 ) 2 + (Y2-Y0) 2} 1/2, is Lp = {(Xp9-X0) 2 + (Yp9-Y0) 2} 1/2.

Next, t = 0.5 is set in step S38-2-2. In step S38-2-3, the point extended from the measurement center of gravity in the direction of the unit vector
Define x = dx × t + X0 and y = dy × t + Y0. Next, in step S38-2-4, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether x <Xs, y <Ys, ImageM (x, y) = 1. judge. If No here, 0.5 is added to t in step S38-2-5. In the case of Yes, in step S38-2-6, Xp4 = x, Yp4 = y and a measurement edge point are obtained (End).

Further, (minus theta _w of the measurement edge points) at step S38-2-11, - calculating the unit vector in the theta _w degree direction. In this embodiment, dx and dy are unit vectors in the −45 degree direction.
dx = {(X2-X0) / L + (Xm9-X0) / Lm} / 2
dy = {(Y2-Y0) / L + (Ym9-Y0) / Lm} / 2
However, L = {(X2-X0 ) 2 + (Y2-Y0) 2} 1/2, is Lm = {(Xm9-X0) 2 + (Ym9-Y0) 2} 1/2.
Next, t = 0.5 is set in step S38-2-12. In step S38-2-13, the point extended from the measurement center of gravity in the direction of the unit vector
Define x = dx × t + X0 and y = dy × t + Y0. Next, in step S38-2-14, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether x <Xs, y <Ys, ImageM (x, y) = 1. judge. If No here, 0.5 is added to t in step S38-2-15. In the case of Yes, in step S38-2-16, Xmp4 = x, Ym4 = y and a measurement edge point are obtained (End).
As described above, the measurement edge points of θ _w degrees and −θ _w degrees can be added.

FIG. 39 is a flowchart showing the flow of the additional processing of the measurement edge point (± second angle) of θ _n degrees and −θ _n degrees (step S38-3, FIGS. 19A and 19D). )reference). Here, an example in which the second angle is 22.5 degrees is shown, but the second angle is an angle narrower than the first angle (10 degrees to 40 degrees).

In step S38-3-1, a unit vector in the +22.5 degree direction is defined as follows.
dx = {(X2-X0) / L + (Xp4-X0) / Lp} / 2
dy = {(Y2-Y0) / L + (Yp4-Y0) / Lp} / 2
However, L = {(X2-X0 ) 2 + (Y2-Y0) 2} 1/2, is Lm = {(Xm4-X0) 2 + (Ym4-Y0) 2} 1/2.

Next, t = 0.5 is set in step S38-3-2. In step S38-3-3, the point extended from the measurement center of gravity in the direction of the unit vector
Define x = dx × t + X0 and y = dy × t + Y0. Next, in step S38-3-4, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether x <Xs, y <Ys, ImageM (x, y) = 1. judge. If No here, 0.5 is added to t in step S38-3-5. In the case of Yes, in step S38-3-6, Xp2 = x, Yp2 = y and an edge point are obtained (End).

In step S38-3-11, a unit vector in the -22.5 degree direction is defined as follows.
dx = {(X2-X0) / L + (Xm4-X0) / Lm} / 2
dy = {(Y2-Y0) / L + (Ym4-Y0) / Lm} / 2
However, L = {(X2-X0 ) 2 + (Y2-Y0) 2} 1/2, is Lm = {(Xm4-X0) 2 + (Ym4-Y0) 2} 1/2.

Next, in step S38-3-12, t = 0.5 is set. In step S38-3-13, the point extended from the measurement center of gravity in the direction of the unit vector
Define x = dx × t + X0 and y = dy × t + Y0. Next, in step S38-3-14, it is determined whether or not the target pixel (x, y) is on the contour, that is, whether x <Xs, y <Ys, ImageM (x, y) = 1. judge. If No here, in step S38-3-15, 0.5 is added to t. In the case of Yes, in step S38-3-16, Xm2 = x, Ym2 = y and a measurement edge point are obtained (End).

Next, a method for calculating the arc area at the rear pole will be described. Area ratio, the calculation of a wide area ranging from - [theta] _w of the + theta _w degree is carried out in two stages in the calculation of the narrow area in the range of - [theta] _n of the + theta _n degrees, first, wide A calculation method in the area will be described.

  FIG. 40A is a diagram in the case of a left direction projection image, and FIG. 40B is a table showing a method for calculating a ratio based on the obtained area values in a wide region.

FIG. 41 is a flowchart showing an example of calculation processing of the arc area in the rear pole portion described with reference to FIG. First, in step S40c-1, the upper area Ap = 0 and v = 0 are initialized, and in step S40c-2, u = 0 is initialized. Next, in step S40c-3,
x = (Xp4-X0) × u / Lup + (X2-X0) × v / Lv
y = (Yp4-Y0) × u / Lup + (Y2-Y0) × v / Lv
X and y in the region surrounded by the +45 degree direction and the 0 degree direction are obtained by the following equation.

０でない場合には、同一画素が既に重複してカウントされていないかをBuff(x,y):0で確認し、Buff(x,y)＝0の場合のみBuff(x,y)=1、Ap¬Ap+1を代入して、ステップＳ４０ｃ−７に進み、u¬u+0.5とする。ステップＳ４０ｃ−８で、uとLupの大小を比較し、u<Lupであれば、ステップＳ４０ｃ−３に戻り、u³Lupであれば、ステップＳ４０ｃ−９に進む。ｖについて判定をし、その結果、u³Lupであれば、上部面積Ａｐが決定される。同様に処理によりＳ４０ｃ−１１以降の処理でＡｍについて求める。Ａｍの算出において、Ａｐの算出と異なる点は、Ｓ４０ｃ−１５、Ｓ４０ｃ−１６の処理で、Ａｐの算出において同一画素が既に重複してカウントされている場合は、ＡｐとＡｍの双方に重複する境界領域であることを示すため、Ａｍ側で二重カウントを抑止するだけでなく、Ｓ４０ｃ−１６の処理で、Ａｐ側に既にカウントされた当該１画素を取り消す処理を行う。このようにして、図２０に示す後極部における円弧面積を求めることができる。 In step S40c-4, it is determined whether (x, y) is within the area calculation target region, that is, ImageM (x, y): 0.

If it is not 0, Buff (x, y): 0 is used to check whether the same pixel has already been counted twice. Only when Buff (x, y) = 0, Buff (x, y) = 1 , Ap¬Ap + 1 is substituted, and the process proceeds to step S40c-7 to set u¬u + 0.5. Step S40C-8, and compares the u and Lup, if u <Lup, the process returns to step S40C-3, if u ³ Lup, the process proceeds to step S40C-9. If v is determined and u ³ Lup is determined as a result, the upper area Ap is determined. Similarly, Am is obtained in the processing after S40c-11 by the processing. In the calculation of Am, the difference from the calculation of Ap is that in the processes of S40c-15 and S40c-16, if the same pixel has already been counted repeatedly in the calculation of Ap, both Ap and Am overlap. In order to indicate the boundary region, not only the double counting is suppressed on the Am side, but also the process of canceling the one pixel already counted on the Ap side is performed in the process of S40c-16. In this manner, the arc area at the rear pole portion shown in FIG. 20 can be obtained.

  Next, a calculation method in a narrow area will be described. Fig.42 (a) is a figure which shows the calculation method (narrow area | region) of the circular arc area in a back pole part. This figure is an example in the case of a left direction projection image. Calculation is performed using Xp2, Yp2 and Xm2, Ym2 obtained above. As shown in FIG. 42B, the E / N ratio and the T / B ratio in the left direction image, the right direction image, the upward direction image, and the downward direction image can be obtained.

  FIG. 43 is a flowchart showing a method of calculating the arc area in the rear pole part (narrow region). It can be calculated in the same manner as the wide area described above. The unit vector is defined above (Formula 1).

Here, the unit vector Up → = (Xp2-X0, Yp2-0Y) / Lup, Um → = (Xm2-X0, Ym2-) with respect to the +22.5 degree direction, -22.5 degree direction, and 0 degree direction. 0Y) / Lum, V → = (X2-X0, Y2-Y0) / Lv.
However, Lup = {(Xp2-X0) ² + (Yp2-Y0) ² } ^1/2 , Lum = {(Xm2-X0) ² + (Ym2-Y0) ² } ^1/2 , Lv = {(X2- X0) ² + (Y2-Y0) ² } ^1/2

  Projected image data after contour extraction is ImageM (x, y) = {0: Background, 1: Contour, 2: Inside} (x = 0,…, Xs; y = 0,…, Ys) Xs = Ys = 300 [pixel] and an image buffer Buff (x, y) of the same size are defined, and the initial values of the image buffer are all set to 0. The image buffer Buff (x, y) prevents the same pixel from being counted twice.

  The upper area is Ap2, the lower area is Am2, and the initial values are 0 respectively. The area ratio is calculated using the calculated Ap2 and Am2.

First, Ap2 = 0 and v = 0 are set in step S39-1, and u = 0 is set in step S39-2. v and u are parameters. In step S39-3, calculation of the left and right upper circular arc areas in the rear pole part of FIG. X and y in the region surrounded by the +22.5 degree direction and the 0 degree direction are obtained.
x = (Xp2-X0) × u / Lup + (X2-X0) × v / Lv
y = (Yp2-Y0) × u / Lup + (Y2-Y0) × v / Lv

In step S39-4, it is calculated whether (x, y) is within the area calculation target region, that is, whether it is ImageM (x, y): 0. In the case of No, in step S39-5, it is confirmed by Buff (x, y): 0 whether or not the same pixel has already been counted, and if it is 0, the same value is stored in the buffer in step S39-6. Buff (x, y) = 1 is stored, and 1 is added to the area Ap2 on the upper side in the horizontal direction. Then, in the case where ImageM (x, y) = 0 in step S39-4, 0.5 is added to u in step S39-8, and u and Lup are compared in step S39-8. Here, if u <Lup, the flow returns to step S39-2, if u ³ Lup, the process proceeds to step S39-9, adding 0.5 to v. In step S39-10, v and Lv are compared. If v <Lv, the process returns to step S39-2. If v ³ Lv, Am2 = 0 and v = 0 in step S39-11. And u = 0 in step S39-12. Thus, Ap2 is obtained.

Next, in step S39-13, x and y in the region surrounded by the -22.5 degree direction and the 0 degree direction are obtained.
x = (Xm2-X0) × u / Lum + (X2-X0) × v / Lv
y = (Ym2-Y0) × u / Lum + (Y2-Y0) × v / Lv
Thus, the arc areas on the lower left and right sides in the rear pole portion of FIG.

Next, in step S39-14, it is determined whether or not (x, y) is within the area calculation target region, that is, whether or not ImageM (x, y) is 0. ImageM (x, y)> If 0, the process proceeds to step S39-15, and it is confirmed by Buff (x, y): 0 whether the same pixel has already been counted twice. If Buff (x, y) is 0, go to Step S39-17,
Buff (x, y) = 2
Am2 ← Am2 + 1
age,
If Buff (x, y) is not 0, go to Step S39-16,
Buff (x, y) = 2
Ap2 ← Ap2-1
And This removes the count of the boundary portion already counted on the upper side.

In step S39-18, 0.5 is added to u, and u and Lum are compared in step S39-19. If u <Lum, the process returns to step S39-12. If u ³ Lum, 0.5 is added to v in step S39-20. In step S39-21, v and Lv are compared. If v <Lv, the process returns to step S39-12. If v ³ Lv, the process is terminated. Thus, Am2 is obtained. As a result, it is possible to calculate by counting the pixels so that the arc area at the rear pole portion of FIG. 20 does not overlap as Ap2 and Am2.

  FIG. 44 is a diagram illustrating an example of a method of calculating an angle formed by the rear pole 3 edge points (sharpness).

Note that, unlike the calculation of the area ratio described above, the calculation of the angle is performed only in the region corresponding to the area of the narrow region.
The analysis processing algorithm has been described above.

Next, the actually obtained projection image and measurement results will be described.
FIG. 45 is an example of a six-direction projection image to be analyzed, and it can be seen that projection images from six directions, ie, left, right, up, down, and front and rear, are obtained. FIG. 46 is a diagram illustrating an example of an analysis parameter setting screen in analysis. As setting parameters, a binarization threshold (range 0 to 255, normal setting value: 25), a deleted spot area (unit: pixel, normal setting value: 200), search range weight (front pole, rear pole) Normal setting values: 16 and 4), lens edge measurement range (range in which Pc is moved in FIGS. 15A (2) and (3), normal setting value: 2 to 12), centroid search range (FIG. 15A (3 ) In which Pg is moved, unit%, normal set value: 30), measurement origin (= measurement center of gravity, distance from back pole edge point, normal set value: 87 pixels), symmetry calculation angle (Area calculation angle range of wide area, normal setting value: 45 degrees), Sharpness calculation angle (sharpness calculation and area calculation angle range of narrow area, normal setting value: 22 degrees), analysis eyeball is left eye Standard eye length upper limit (unit: pixel, common eye) Setting value: 175 pixels), left-right symmetry lower limit value (%, normal setting value: 90%), upper limit value (%, normal setting value: 110%), vertical symmetry lower limit value (%, normal setting) Value: 90%), upper limit (%, normal setting value: 110%), rear pole sharpness lower limit value (normal setting value: 140 degrees), upper limit value (normal setting value: 150 degrees) In general, however, it is only necessary to set whether the analysis eyeball is the left eye or the right eye.

  FIG. 47 is a diagram illustrating an example of the analysis result. Analysis results based on the above analysis process can be viewed on a single screen. For example, for L (left projection image), R (right projection image), T (upward projection image), and B (downward projection image), the axial length (unit: pixel), rear polar sharpness (0 ~ 360 degrees), narrow area ratio (%), wide area ratio (%), area ratio (T / B, E / N, one of the values of narrow area or wide area area ratio, difference from 100 Large value), rear pole peak position (coordinate value, pixel unit), left-right symmetry (ear, left-right symmetry, nose-side judgment result), vertical symmetry (upward, left-right symmetry, bottom) The determination result of any one of the shifts, the rear sharp tip type (the determination result of any of barrel type, intermediate type, spindle type), and the type classification (four digits in Table 1) are displayed. The type classification is 2000, the rear pole sharp type is barrel type 2, and the other axial length, left-right direction symmetry, and up-down direction symmetry are normal. From this screen, the analysis results and pathological conditions can be seen at a glance.

  FIG. 48 is a diagram corresponding to FIG. 45 and is a diagram illustrating an example of a confirmation screen for the auxiliary line setting status at the time of analysis. FIG. 48 shows auxiliary lines (straight lines connecting various measurement edge points from the measurement center of gravity) for correcting the eye axis in the projected images of the upper, lower, left, and right eyeballs. On this confirmation screen, the eye axis can be corrected and reanalysis can be performed with the designated auxiliary line.

The above is an example in the normal case, but an example in the case of intense myopia is shown below.
FIG. 49 is a diagram illustrating an example of a six-direction projection image to be analyzed in the case of intensity myopia. It can be seen that the image is quite different from the image shown in FIG. FIG. 50 shows an example of an analysis parameter setting screen, which is the same as FIG. An example of the result of analysis with the same setting will be described with reference to FIG. It can be seen that the analysis result is quite different from FIG. In particular, the axial length, rear pole sharpness, left-right symmetry, vertical symmetry, rear pole sharpness, and the like are different. In the type classification, it can be seen in Table 1 that the type classification is 11, the axial length is abnormal, and the left-right symmetry is near the nose (less than 90%). FIG. 52 shows a confirmation screen for the auxiliary line setting status at the time of analysis. FIG. 52 shows auxiliary lines for correcting the eye axis in the projected images of the eyeballs before and after the upper, lower, left, and right. On this confirmation screen, the eye axis can be corrected and reanalysis can be performed with the designated auxiliary line.

  FIG. 53 is a diagram illustrating an example of a 6-direction projection image of a 6-direction projection image (mild) to be analyzed. Compared to FIG. 49, it can be seen that the projected image is somewhat similar to FIG. In this case, as shown in FIG. 54, the left eye is the analysis target. The result is FIG. 55 and the type classification is 121. The axial length is abnormal, the left-right direction symmetry is closer to the ear (110% or more), and the vertical direction symmetry is lower (less than 90%) )It can be seen that it is.

  FIG. 56 is a diagram showing an example of a confirmation screen for the auxiliary line setting status at the time of analysis. As shown in FIG. 56, regarding the automatic setting of the eye axis, it is clear that it is out of the parallel / vertical direction of the screen and is clearly inappropriate. FIG. 57 is a confirmation screen for the auxiliary line setting status at the time of analysis, and shows a state in which correction of the eye axis (manual correction in this example) is added. As shown in FIG. 57, the auxiliary line of the eye axis coincides with the parallel / vertical direction of the screen, and it can be confirmed that it is corrected correctly.

  FIG. 58 is an analysis parameter setting screen. For example, the setting value is not changed from the setting screen shown in FIG. 54, but the setting value can be confirmed on the screen.

  FIG. 59 is a diagram illustrating an example of a reanalysis result after correction. Compared with FIG. 55, the symmetry in the left-right direction has changed from the near side to the left-right symmetry, and the posterior pole sharpness has changed from an intermediate type to a spindle type (normal eyes are barrel-shaped as shown in FIG. 47). ), The type classification has changed from 121 to 1101. Thus, the correct value can be obtained by obtaining the analysis result after changing the eye axis.

  In the above-described embodiment, the configuration and the like illustrated in the accompanying drawings are not limited to these, and can be appropriately changed within a range in which the effect of the present invention is exhibited. In addition, various modifications can be made without departing from the scope of the object of the present invention. It may be a program for causing a computer to execute the image processing method of the eyeball MRI, or a computer-readable recording medium for recording the program.

(Summary)
According to the present invention, it is possible to contribute to an early diagnosis of “intensity myopia” without imposing a burden on the patient, and to provide objective diagnosis data that does not depend on the subjectivity of the interpretation doctor.

  Further, it can be realized by inexpensive general-purpose personal computer software, and it is not necessary to add a special graphics board such as a GPU (Graphics Processing Unit), so that an inexpensive diagnostic service can be provided.

  In addition, without introducing new ophthalmology equipment for the image input system, existing popular MRI equipment can be diverted and radiation is not used, so there is no problem of exposure, and it is possible to make repeated diagnoses regularly. There is no need for customization. That is, according to the present invention, an inexpensive system can be constructed without investing in new hardware equipment.

(Appendix)
The present invention includes the following inventions.
(1) a first step of making a diagnosis with an ophthalmologic diagnostic device;
If diagnosis is not possible in the first step, MRI imaging of the head and extracting a 3D voxel image of the eyeball part;
Rendering the projection image data in any one of the six directions of up / down / left / right / front / back based on the three-dimensional voxel image of the eyeball part;
Based on at least one of the parameters of the eyeball shape by image analysis, the rear axis sharpness, the left / right symmetry, the left / right symmetry, the up / down symmetry, and the number of rear pole protrusions with respect to the projected image. Diagnosing myopia; and
A method for diagnosing intense myopia, comprising:
(2) The method for diagnosing intense myopia according to (1), further comprising a step of performing classification processing of the type of intense myopia for each shape parameter.

  As described above, when the diagnosis of intense myopia is not confirmed in a general ophthalmic diagnosis, it is preferable to use the image analysis technique for ophthalmic diseases using head MRI imaging data according to the present embodiment.

  The present invention can be used in an image analysis apparatus for ophthalmic diseases.

DESCRIPTION OF SYMBOLS A ... Image analysis apparatus, 1 ... Image input part, 3 ... Voxel image extraction part, 5 ... Projection image calculation part, 7 ... Projection image analysis part, 11 ... Output part, 7-1 ... Asymmetry calculation part, 7-2 ... sharpness calculation part, 7-3 ... rear pole protrusion peak calculation part, 7-4 ... eye axis automatic setting part, 7-5 ... vertex correction part.

Claims (17)

An image input unit for acquiring a three-dimensional voxel image of the head including the eyeball;
A voxel image extraction unit that extracts a three-dimensional eyeball voxel image of at least one eye of the eyeball region for each of the left and right eyeballs from the input three-dimensional voxel image;
A projection image calculation unit that creates a first projection image obtained by projecting the extracted three-dimensional eyeball voxel image in at least one of the upper and lower directions;
Setting an eye axis connecting the anterior ocular pole vertex and the posterior ocular pole vertex to the first projection image, and at least the ocular axial length between the ocular anterior pole vertex and the ocular posterior pole vertex; A projection image analysis unit that calculates the degree of asymmetry in the left-right direction based on a region obtained by dividing the posterior pole side into the left and right with the eye axis as a boundary, and an image analysis parameter;
An image analysis apparatus for ophthalmic diseases, comprising:   The ophthalmic disease image analysis apparatus according to claim 1, wherein the three-dimensional voxel image is acquired using MRI. The projection image calculation unit
Creating a second projection image obtained by projecting the three-dimensional eyeball voxel image in at least one of the left and right directions;
The image analysis unit sets an eye axis with respect to the second projection image, and sets the rear pole side at the boundary between the eye axis length and the eye axis between at least the anterior anterior vertex of the eyeball and the posterior apex of the eyeball. The image analysis apparatus for ophthalmic diseases according to claim 1 or 2, wherein the degree of asymmetry in the vertical direction is calculated as a parameter for image analysis based on the vertically divided region. The image analysis unit
The sharpness degree of the arcuate curvature of the rear pole is calculated as an image analysis parameter for at least one of the first and second projection images. The image analysis apparatus for ophthalmic diseases according to claim 1. The image analysis unit
5. The calculation result of at least one of the first and second projection images is a value obtained by averaging at least one of the left and right or upper and lower sets of parameters. The image analysis apparatus for ophthalmic diseases according to any one of the above. The projection image calculation unit
Creating a third projected image projected backward with respect to the three-dimensional eyeball voxel image;
The said image analysis part calculates the number of the peak points corresponding to a back pole protrusion part with respect to a said 3rd projection image as a parameter of image analysis. The image analysis apparatus for ophthalmic diseases according to item 1.   The image analysis unit extracts an outline of an eyeball from the third projection image, further sets a center of gravity of the eyeball, and sets the position of each calculated peak point as an image analysis parameter. The image analysis apparatus for ophthalmic diseases according to claim 6, wherein the image is calculated based on a distance from the center of gravity, which is a polar coordinate based on the center of gravity point, and an angle in a peak point direction from the center of gravity. The image analysis unit
Contour line extraction of the eyeball part is performed on the projected image, and the anterior ocular vertex and the posterior ocular vertex are located on the ocular outline, and the anterior ocular vertex and the posterior eyeball The straight line connecting the extreme vertices is equal in distance to the intersection with the contour line extending vertically from side to side or up and down from several vertices of the straight line, and The image analysis apparatus for ophthalmic diseases according to any one of claims 1 to 3, wherein the eye axis is automatically set by setting a vertex of the posterior pole of the eyeball. The image analysis unit
Extracting the outline of the eyeball from the projected image, further setting the center of gravity of the eyeball, the apex of the front of the eyeball and the apex of the back of the eyeball being located on the outline of the eyeball, and The eyeball front pole vertex, the eyeball rear pole vertex, and the barycentric point are set so that a polygonal line connecting the three points of the eyeball front pole vertex, the center of gravity point, and the eyeball rear pole vertex is as straight as possible. The image analysis apparatus for ophthalmic diseases according to any one of claims 1 to 3, wherein the eye axis is automatically set by setting. The image analysis unit
Contour extraction of the eyeball part with respect to the first or second projection image, and a dialog that is set in an interactive form with the screen for either the frontal vertex of the eyeball or the rearward vertex of the eyeball The contour line on the extension line of the straight line connecting the anterior ocular pole vertex and the posterior pole apex with respect to the anterior ocular pole vertex or the ocular posterior pole apex set by the mold apex input unit The image analysis apparatus for ophthalmic diseases according to any one of claims 1 to 3, further comprising a vertex correction unit that automatically corrects the intersection of the two. The image analysis unit
Contour extraction of the eyeball part is performed on the first or second projection image, a measurement center of gravity point is set in the eyeball part, and a reference line connecting the vertex of the back of the eyeball from the measurement center of gravity point is set. Then, a straight line A and a straight line B in the direction of the + first angle and the first angle are set around the center of gravity for measurement, and the area inside the contour line surrounded by the straight line A and the reference line 2. An area B inside the contour line surrounded by A, the straight line B, and the reference line is calculated, and a ratio of the area A to the area B is given as the degree of asymmetry. The image analysis apparatus for ophthalmic diseases according to any one of items 1 to 3. The image analysis unit
Contour extraction of the eyeball part is performed on the first or second projection image, a measurement center of gravity point is set in the eyeball part, and a reference line connecting the vertex of the back of the eyeball from the measurement center of gravity point is set. Then, a straight line A and a straight line B are set in the direction of the second angle and the second angle around the center of gravity for measurement, and the intersection A between the straight line A and the contour line and the vertex of the posterior pole of the eyeball The image analysis apparatus for ophthalmic diseases according to claim 4, wherein an angle formed by three vertices of the intersection B of the straight line B and the contour line is given as the sharpness degree.   The image analysis unit has a straight line C and a straight line in the direction of the second angle and the second angle with respect to the reference line connecting the center of gravity for measurement to the apex of the posterior pole of the eyeball. D is additionally set, and an area C inside the contour line surrounded by the straight line C and the reference line and an area D inside the contour line surrounded by the straight line D and the reference line are additionally calculated. The image analysis apparatus for ophthalmic diseases according to claim 11, wherein the degree of asymmetry is a larger value of the ratio of the area C to the area D and the ratio of the area A to the area B. The image analysis unit
The center of gravity point for measurement is provided on the eye axis connecting the apex of the front of the eyeball and the apex of the back of the eyeball at a predetermined distance from the apex of the back of the eyeball. The image analysis apparatus for ophthalmic diseases according to any one of 11 to 13. The image analysis unit
The ophthalmic disease image analysis apparatus according to claim 14, wherein the predetermined distance is set to ½ of an average axial length of an eyeball of a normal person. An image input step of acquiring a three-dimensional voxel image of the head including the eyeball;
A voxel image extraction step of extracting a three-dimensional eye voxel image of at least one eye of the eyeball region for each of the right and left eyeballs from the input three-dimensional voxel image;
A projection image calculation step of creating a first projection image obtained by projecting the extracted three-dimensional eyeball voxel image in at least one of the upper and lower directions;
Setting an eye axis connecting the anterior ocular pole vertex and the posterior ocular pole vertex to the first projection image, and at least the ocular axial length between the ocular anterior pole vertex and the ocular posterior pole vertex; A projection image analysis step for calculating the degree of asymmetry in the left-right direction based on a region obtained by dividing the rear pole side into the left and right with the eye axis as a boundary;
An image analysis method for ophthalmic diseases, comprising:   A program for causing a computer to execute the image analysis method for ophthalmic diseases according to claim 16.

Images