Adjusting and Improving X-ray image's contrast and its quality

Question

I want to adjust/improve the contrast of X-ray images of the hand. I have about 10000 of these X-ray images. Instead of manually editing them, I want to code or automate them. Most of the images in the dataset have image qualities like these three images..

I've tried suggestions by here, and here. However, when I run some of image samples through, I get the same output as the input.

Are there any better ways to improve the contrast particularly for hand X-ray images? So, from these three image inputs, I want them to look like these. If can be done automatically then, that would be awesome too. How would one code this?

Check out the techniques [in this page](https://docs.opencv.org/4.x/d5/daf/tutorial_py_histogram_equalization.html) — Jeru Luke, Jul 20 '22 at 20:16
Try CLAHE with different parameters. You may also try Bilateral Filtering before CLAHE (for reducing noise) or [BM3D](https://pypi.org/project/bm3d/). The issue with your question is that you are asking for an algorithm, and it is too general. — Rotem, Jul 20 '22 at 20:27

fmw42 · Accepted Answer · 2022-07-21T23:42:53.157

I would suggest using skimage (rescale_intensity) to stretch the dynamic range (after covering over the label with the average color) in Python/OpenCV/Skimage. It automatically finds the input min and max values and stretches those to full black and full white (i.e. 0 and 255) or compute the min and max values and bias if desired.

Read the input
Convert to grayscale
blur
Apply morphology close
Threshold
Get contours and exclude too small ones and too large ones so that one selects the label.
Draw a filled contour on a copy of the grayscale image of value equal to the mean of the grayscale image in order to cover over the label
Use Skimage to stretch the min and max values to 0 and 255 respectively. If you do not need to bias the min or max, then you can replace (minval,maxval) with 'image'. It will automatically compute the minval and maxval and you will not need to use Numpy to find them.
Save the output

import cv2
import numpy as np
import skimage.exposure

# load images
img1 = cv2.imread('xray1.webp')
img2 = cv2.imread('xray2.webp')
img3 = cv2.imread('xray3.webp')

# convert to gray
gray1 = cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)
gray2 = cv2.cvtColor(img2,cv2.COLOR_BGR2GRAY)
gray3 = cv2.cvtColor(img3,cv2.COLOR_BGR2GRAY)

# blur
blur1 = cv2.GaussianBlur(gray1, (0,0), sigmaX=6, sigmaY=6)
blur2 = cv2.GaussianBlur(gray2, (0,0), sigmaX=6, sigmaY=6)
blur3 = cv2.GaussianBlur(gray3, (0,0), sigmaX=6, sigmaY=6)

# morphology
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (45,45))
morph1 = cv2.morphologyEx(blur1, cv2.MORPH_CLOSE, kernel)
morph2 = cv2.morphologyEx(blur2, cv2.MORPH_CLOSE, kernel)
morph3 = cv2.morphologyEx(blur3, cv2.MORPH_CLOSE, kernel)

# threshold
thresh1 = cv2.threshold(morph1, 0, 255, cv2.THRESH_OTSU)[1] 
thresh2 = cv2.threshold(morph2, 0, 255, cv2.THRESH_OTSU)[1] 
thresh3 = cv2.threshold(morph3, 0, 255, cv2.THRESH_OTSU)[1] 

# get contours and filter on size
masked1 = gray1.copy()
meanval = int(np.mean(masked1))
contours = cv2.findContours(thresh1, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contours = contours[0] if len(contours) == 2 else contours[1]
for cntr in contours:
    area = cv2.contourArea(cntr)
    if area > 500 and area < 50000:
        cv2.drawContours(masked1, [cntr], 0, (meanval), -1)
    
masked2 = gray2.copy() 
meanval = int(np.mean(masked2))
contours = cv2.findContours(thresh2, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contours = contours[0] if len(contours) == 2 else contours[1]
for cntr in contours:
    area = cv2.contourArea(cntr)
    if area > 500 and area < 50000:
        cv2.drawContours(masked2, [cntr], 0, (meanval), -1)

masked3 = gray3.copy() 
meanval = int(np.mean(masked3))
contours = cv2.findContours(thresh3, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contours = contours[0] if len(contours) == 2 else contours[1]
for cntr in contours:
    area = cv2.contourArea(cntr)
    if area > 500 and area < 50000:
        cv2.drawContours(masked3, [cntr], 0, (meanval), -1)

# stretch
minval = int(np.amin(masked1))
maxval = int(np.amax(masked1))
result1 = skimage.exposure.rescale_intensity(masked1, in_range=(minval,maxval), out_range=(0,255)).astype(np.uint8)
minval = int(np.amin(masked2))
maxval = int(np.amax(masked2))
result2 = skimage.exposure.rescale_intensity(masked2, in_range=(minval,maxval), out_range=(0,255)).astype(np.uint8)
minval = int(np.amin(masked3))
maxval = int(np.amax(masked3))
result3 = skimage.exposure.rescale_intensity(masked3, in_range=(minval,maxval), out_range=(0,255)).astype(np.uint8)

# save output
cv2.imwrite('xray1_stretched.png', result1)
cv2.imwrite('xray2_stretched.png', result2)
cv2.imwrite('xray3_stretched.png', result3)

# Display various images to see the steps
cv2.imshow('thresh1', thresh1)
cv2.imshow('thresh2', thresh2)
cv2.imshow('thresh3', thresh3)
cv2.imshow('result1', result1)
cv2.imshow('result2', result2)
cv2.imshow('result3', result3)
cv2.waitKey(0)
cv2.destroyAllWindows()

Result 1:

Result 2:

Result 3:

I was wondering, have you tried global or local equalization? — Jeru Luke, Jul 21 '22 at 14:42
No, I have not tried using CLAHE if that is what you are asking on these images. But I was concerned about using them if they clip data. This process that I used limits the stretch to full dynamic range and no more unless one wants to bias the minval or maxval. Feel free to post examples of equalization. — fmw42, Jul 21 '22 at 15:32
Thank you for the reply and answer @fmw42. When you mentioned: "If you do not need to bias the min or max, then you can replace (minval,maxval) with 'image'", do you mean in the code where you mention `in_range=(minval,maxval)`, it should become like `in_range = img1`? — Bathtub, Jul 23 '22 at 15:01
`result1 = skimage.exposure.rescale_intensity(masked1, in_range='image', out_range=(0,255)).astype(np.uint8)` See docs at https://scikit-image.org/docs/stable/api/skimage.exposure.html#skimage.exposure.rescale_intensity — fmw42, Jul 23 '22 at 16:19

score 1 · Answer 2 · answered Jul 22 '22 at 06:58

I tried a little to a local-type-enhancement (with C++).
Looking results, I'm worried about strong halo...

Note : I tested with 50% size images, because images are too large for my monitor.

class MyTryImpl
{
public:
    MyTryImpl( double Sigmoid_A=5, unsigned char T_min=5, unsigned char T_max=250 )
    {   MakeTable( Sigmoid_A, T_min, T_max );   }

public:
    void Proc( const cv::Mat &Src8UC1, const cv::Mat &GlobalLMap, cv::Mat &Dst8UC1 ) const
    {
        if( Src8UC1.empty() || Src8UC1.type()!=CV_8UC1 ){   throw std::invalid_argument( "must : Src is 8UC1" );    }
        if( GlobalLMap.empty() || GlobalLMap.type()!=CV_8UC1 ){ throw std::invalid_argument( "must : GlobalLMap is 8UC1" ); }
        if( Src8UC1.size() != GlobalLMap.size() ){  throw std::invalid_argument( "must : Src8UC1.size() == GlobalLMap.size()" );    }

        Dst8UC1.create( Src8UC1.size(), CV_8UC1 );
        for( int y=0; y<Src8UC1.rows; ++y )
        {
            const unsigned char *pS = Src8UC1.ptr<unsigned char>(y);
            const unsigned char *pB = GlobalLMap.ptr<unsigned char>(y);
            unsigned char *pR = Dst8UC1.ptr<unsigned char>(y);
            for( int x=0; x<Src8UC1.cols; ++x, ++pS,++pB,++pR )
            {   *pR = m_ResultTbl[ *pS ][ *pB ];    }
        }
    }

    void Proc(
        const cv::Mat &Src8UC1, cv::Mat &Dst8UC1,
        int BlurKernelSize,
        bool WithGaussian=true  //if ture GaussianFilter, else BoxFilter
    ) const
    {
        if( Src8UC1.empty() || Src8UC1.type()!=CV_8UC1 )
        {   throw std::invalid_argument( "must : Src is 8UC1" );    }

        cv::Mat Blurred;
        {
            int s = std::max( 3, BlurKernelSize | 0x01 );
            if( WithGaussian ){ cv::GaussianBlur( Src8UC1, Blurred, cv::Size(s,s), 0 ); }
            else {  cv::blur( Src8UC1, Blurred, cv::Size(s,s) );    }
        }
        Proc( Src8UC1, Blurred, Dst8UC1 );
    }

private:
    void MakeTable( double Sigmoid_A, unsigned char T_min, unsigned char T_max )
    {
        if( T_min==0 ){ throw std::invalid_argument( "must : T_min > 0" );  }
        if( T_max==255 ){   throw std::invalid_argument( "must : T_max < 255" );    }
        if( T_min>=T_max ){ throw std::invalid_argument( "must : T_min < T_max" );  }

        unsigned char B = 0;
        while( true )
        {
            double b = std::max( T_min, std::min(T_max,B) ) / 255.0;
            const double g = log(0.5) / log( b<=0.5 ? b : 1.0-b);
            unsigned char S=0;
            while( true )
            {
                double s = S / 255.0;
                double c = Sig( Gam( (b<=0.5 ? s : 1.0-s), g ), Sigmoid_A );
                m_ResultTbl[S][B] = cvRound( 255.0 * (b<=0.5 ? c : 1.0-c) );
                if( S==0xFF )break;
                ++S;
            }
            if( B==0xFF )break;
            ++B;
        }
    }

    static double Sig( double x, double a )
    {
        double exp1 = exp( -a*(2*x -1) );
        double exp2 = exp( -a );
        double nume = (1-exp1)*(1+exp2);
        double denom = (1+exp1)*(1-exp2);
        return 0.5 * ( 1 + nume/denom );
    }

    static double Gam( double x, double g ){    return pow( x, g ); }

private:
    unsigned char m_ResultTbl[256][256];
};

int main()
{//Test for 3 imgs
    const std::string SrcImgFileNames[3] = { "Xray1.png", "Xray2.png", "Xray3.png" };
    const std::string SaveFileNames[3] = { "Result1.png", "Result2.png", "Result3.png" };
    MyTryImpl MyTest;
    for( int i=0; i<3; ++i )
    {
        cv::Mat SrcImg = cv::imread( SrcImgFileNames[i], cv::IMREAD_GRAYSCALE );
        if( SrcImg.empty() )return 0;

        cv::Mat ResultImg;
        MyTest.Proc( SrcImg, ResultImg, SrcImg.cols/10, true );
        cv::imwrite( SaveFileNames[i], ResultImg );
    }
    return 0;
}

Results:

With this method, local average value in source image becomes to middle level(about 128) in result image, so background became to gray. — fana, Jul 22 '22 at 07:06

Adjusting and Improving X-ray image's contrast and its quality

2 Answers2