Image editing using ArrayFire: Part 3

Today, we will be doing the third post in our series Image editing using ArrayFire. References to old posts are available below.
* Part 1
* Part 2

In this post, we will be looking at the following operations.

• Image Histogram
• Simple Binary Theshold
• Otsu Threshold
• Iterative Threshold
• Adaptive Binary Threshold
• Emboss Filter

Today’s post will be mostly dominated by different types of threshold operations we can achieve using ArrayFire.

Image Histogram

We have a built-in function in ArrayFire that creates a histogram. The input image was converted to gray scale before histogram calculation as our histogram implementation works for vector and 2D matrices only.

Input

In case, you need histogram for all three channels of a color image, you can do it using the following four line code stub.

array color = loadimage("",true); // last parameter is to indicate that you would like 
                                             // array to have all channels of the image
array hist(256,color.dims(2)); // 2nd dimension is our channel dimension
                               // ArrayFire stores channels as independent slices currently.
hist(span,0) = histogram(color(span,span,0),256,0,255);
hist(span,1) = histogram(color(span,span,1),256,0,255);
hist(span,2) = histogram(color(span,span,2),256,0,255);

Simple Binary Threshold

A simple threshold operation can be achieved very easily in ArrayFire using the awesome element wise operations. For this particular use case, we will use the greater than and less than operators.

array threshold(const array &in, float thresholdValue)
{
    int channels = in.dims(2);
    array ret_val= in.copy();
    if (channels>1)
        ret_val = colorspace(in,af_gray,af_rgb);
    ret_val = 255.0f*(ret_val>thresholdValue);
    return ret_val;
}

Let’s go through what’s going on in this function line by line. We initially check if the image has more than one channel and convert it to gray scale image if required. Following that, we set the pixels that are less than threshold to 0 and the rest of them to 255 255.0f*(ret_val>thresholdValue). If you are interested in knowing about what other element wise operations you can do on the objects of type array, you can find the details here.

Input

Simple binary threshold

Otsu Threshold

The famous Otsu threshold method is very reliable in the case of images with bi-modal histograms. Otsu method calculates a global threshold value using the histogram of the image rather than working on the image itself. Since the histogram is only a array of 256 elements, we chose to do the threshold computation on CPU. ArrayFire is used to calculate the histogram and finally threshold the image using the threshold value computed by Otsu method. The image threshold operation reuses the simple threshold function discussed in the previous section. Below given is a comparison of simple binary threshold and Otsu threshold for your reference.

/**
 * Note:
 * suffix B indicates subset of all graylevels before current gray level
 * suffix F indicates subset of all gray levels after current gray level
 */
array otsu(const array &in)
{
    array gray;
    int channels = in.dims(2);
    if (channels>1)
        gray  = colorspace(in,af_gray,af_rgb);
    else
        gray  = in;
    unsigned total = gray.elements();
    array hist  = histogram(gray,256,0.0f,255.0f);
    array wts   = array(seq(256));

    array wtB   = accum(hist);
    array wtF   = total - wtB;
    array sumB  = accum(wts*hist);
    array meanB = sumB/wtB;
    array meanF = (sumB(255)-sumB) / wtF;
    array mDiff = meanB-meanF;

    array interClsVar= wtB * wtF * mDiff * mDiff;
    float max        = af::max(interClsVar);
    float threshold2 = where(interClsVar==max).scalar();
    array threshIdx  = where(interClsVar>=max);
    float threshold1 = threshIdx.elements()>0 ? threshIdx.scalar() : 0.0f;

    return threshold(gray,(threshold1+threshold2)/2.0f);
}

Input image with Bimodal histogram

Simple binary threshold

Otsu threshold

Iterative Threshold

This is a global threshold method where the input image is repetitively thresholded. In each iteration, threshold value is calculated from the previous step’s result. In the first iteration, the threshold value is selected randomly. When the difference between two successive steps threshold is less than certain value, the algorithm exits the iterative loop and the last iteration’s threshold value is used to threshold the original input image.

array iterativeThreshold(const array &in)
{
    array ret_val = colorspace(in,af_gray,af_rgb);
    float T = mean(ret_val), prevT=0.0f;
    do {
        array region1 = (ret_val > T)*ret_val;
        array region2 = (ret_val <= T)*ret_val;
        float r1_avg  = mean(region1);
        float r2_avg  = mean(region2);
        prevT = T;
        T = (r1_avg+r2_avg)/2.0f;        
    } while (abs(T-prevT)>0.01f);
    return threshold(ret_val,T);
}

Input

Otsu threhsold

Iterative threshold

Adaptive Binary Threshold

Threshold methods are usually categorized into two types: global and variable. Global threshold methods have a single threshold value that helps classify the image pixels into two classes. Variable threshold is where local neighborhood properties of the current pixel are used to threshold that very pixel. Adaptive binary threshold is where every pixel is thresholded using a value that is computed using the pixel intensities from the spatial neighborhood of the current pixel. In the below implementation, we used the local operators such as Mean, Median and Minimum & Maximum to calculate the pixel threshold from local neighborhood. Threshold output for each local operator mentioned above are provided below for your reference. The algorithm essentially follows theses main steps

1. Convert image to gray scale if required
2. Convolve the gray scale image using the local operator
3. Take a difference between the convolved image and gray scale image
4. Threshold the difference image using constant
5. Invert the thresholded image and return the result

typedef enum {
    MEAN=0,
    MEDIAN,
    MINMAX_AVG
} LocalThresholdType;

array adaptiveThreshold(const array &in, LocalThresholdType kind, int window_size, int constnt)
{
    int wr = window_size;
    array ret_val = colorspace(in,af_gray,af_rgb);
    if (kind==MEAN) {
        array wind = constant(1,wr,wr)/(wr*wr);
        array mean = convolve(ret_val,wind);
        array diff = mean - ret_val;
        ret_val    = 255.f*(diff>constnt);
    } else if (kind==MEDIAN) {
        array medf = medfilt(ret_val,wr,wr);
        array diff = medf - ret_val;
        ret_val    = 255.f*(diff>constnt);
    } else if (kind==MINMAX_AVG) {
        array minf = minfilt(ret_val,wr,wr);
        array maxf = maxfilt(ret_val,wr,wr);
        array mean = (minf+maxf)/2.0f;
        array diff = mean - ret_val;
        ret_val    = 255.f*(diff>constnt);
    }
    ret_val = 255.f - ret_val;
    return ret_val;
}

Input

Adaptive Threshold with Median as local operator

Adaptive Threshold with Minimum & Maximum as local operator

Adaptive Threshold with Mean as local operator

Emboss Filter

The emboss function is tailored after the emboss filter provided in GIMP, an open source image editing software. The final outcome of emboss function is based on its three input parameters: azimuth, elevation & depth. The code and the sample outputs of emboss function are provided below.

/**
 * azimuth range is [0-360]
 * elevation range is [0-180]
 * depth range is [1-100]
 * Note: this function has been tailored after
 * the emboss implementation in GIMP editor
 **/
array emboss(const array &input, float azimuth, float elevation, float depth)
{
    if (depth>100 || depth<1) {
        printf("Depth should be in the range of 1-100");
        return input;
    }
    static float x[3] = {-1,0,1};
    static array hg(3,x);
    static array vg = hg.T();

    array in = input;
    if (in.dims(2)>1)
        in = colorspace(input,af_gray,af_rgb);
    else
        in = input;

    // convert angles to radians
    float phi   = elevation*af::Pi/180.0f;
    float theta = azimuth*af::Pi/180.0f;

    // compute light pos in cartesian coordinates
    // and scale with maximum intensity
    // phi will effect the amount of we intend to put
    // on a pixel
    float pos[3];
    pos[0] = 255.99f * cos(phi)*cos(theta);
    pos[1] = 255.99f * cos(phi)*sin(theta);
    pos[2] = 255.99f * sin(phi);

    // compute gradient vector
    array gx = filter(in,vg);
    array gy = filter(in,hg);

    // create a array from depth value which 
    // is of same dimensions as gradient vector array
    float pxlz  = (6*255.0f)/depth;
    array zdepth= constant(pxlz,gx.dims());
    // dot product of light position vector and gradient vector
    array vdot  = gx*pos[0] + gy*pos[1] + pxlz*pos[2];
    // divide dot product by gradient vector magnitude to
    // calculate scalar projection value
    array norm  = vdot/sqrt(gx*gx+gy*gy+zdepth*zdepth);
    // compute shading condition
    array outwd = vdot < 0.0f;
    // choose shading color based on the above condition
    array color = (1-outwd) * norm;
    return color;
}

Emboss with low elevation angle

Emboss with high elevation value

Emboss with low azimuth angle

Emboss with high azimuth angle

Emboss with low depth value

Emboss with high depth value

Conclusion

All the functions from this post are updated and available through github repository located here.