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?

Answer 1

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:

Answer 2

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: