fast threshold and bit packing algorithm ( possible improvements ? )

Question

I am working on an algorithm that performs a global thresholding of an 8-bit grayscale image into a 1-bit ( bit packed, such that 1 byte contains 8 pixels ) monochrome image. Each pixel in the Grayscale image can have a luminance value of 0 - 255.

My environment is Win32 in Microsoft Visual Studio C++.

I am interested in optimizing the algorithm as much as possible out curiosity, The 1-bit image will be turned into a TIFF. Currently I am setting the FillOrder to be MSB2LSB (Most Significant Bit to Least Significant Bit ) just because the TIFF spec suggests this ( it doesn't necessarily need to be MSB2LSB )

Just some background for those who don't know:

MSB2LSB orders pixels from left to right in a byte just as pixels are oriented in an image as the X coordinate increases. If you are traversing the Grayscale image from left to right on the X axis this obviously requires you to think "backward" as you are packing the bits in your current byte. With that said, let me show you what I currently have (This is in C, I have not attempted ASM or Compiler Intrinsics yet just because I have very little experience with it, but that would be a possibility ).

Because the Monochrome image will have 8 pixels per byte, the Width of the monochrome image will be

(grayscaleWidth+7)/8;

FYI, I assume my largest image to be 6000 pixels wide:

First thing I do (before any image is processed) is

1) calculate a look up table of amounts I need to shift into a specific byte given an X coordinate from my grayscale image:

int _shift_lut[6000];

for( int x = 0 ; x < 6000; x++)
{ 
    _shift_lut[x] = 7-(x%8);
}

With this lookup table I can pack a monochrome bit value into the current byte I am working on with something like:

monochrome_pixel |= 1 << _shift_lut[ grayX ];

which ends up being a huge speed increase over doing

monochrome_pixel |= 1 << _shift_lut[ 7-(x%8)];

The second lookup table I calculate is a Lookup Table that tells me the X index into my monochrome pixel given an X pixel on the Grayscale pixel. This very simple LUT is calculated like such:

int xOffsetLut[6000];
int element_size=8; //8 bits
for( int x = 0; x < 6000; x++)
{
    xOffsetLut[x]=x/element_size;
}

This LUT allows me to do things like

monochrome_image[ xOffsetLut[ GrayX ] ] = packed_byte; //packed byte contains 8 pixels

My Grayscale image is a simple unsigned char* and so is my monochrome image;

This is how I initialize the monochrome image:

int bitPackedScanlineStride = (grayscaleWidth+7)/8;
int bitpackedLength=bitPackedScanlineStride * grayscaleHeight;
unsigned char * bitpack_image = new unsigned char[bitpackedLength];
memset(bitpack_image,0,bitpackedLength);

Then I call my binarize function like such:

binarize(
    gray_image.DataPtr(),
    bitpack_image,
    globalFormThreshold,
    grayscaleWidth,
    grayscaleHeight,
    bitPackedScanlineStride,
    bitpackedLength,
    _shift_lut,  
    xOffsetLut);

And here is my Binarize function ( as you can see I did some loop unrolling, which may or may not help ).

void binarize( unsigned char grayImage[], unsigned char bitPackImage[], int threshold, int grayscaleWidth, int grayscaleHeight, int  bitPackedScanlineStride, int bitpackedLength,  int shiftLUT[], int xOffsetLUT[] )
{
    int yoff;
    int byoff;
    unsigned char bitpackPel=0;
    unsigned char pel1=0;
    unsigned char  pel2=0;
    unsigned char  pel3=0;
    unsigned char  pel4=0;
    unsigned char  pel5=0;
    unsigned char  pel6=0;
    unsigned char  pel7=0;
    unsigned char  pel8=0;
    int checkX=grayscaleWidth;
    int checkY=grayscaleHeight;

    for ( int by = 0 ; by < checkY; by++)
    {
    yoff=by*grayscaleWidth;
    byoff=by*bitPackedScanlineStride;

    for( int bx = 0; bx < checkX; bx+=32)
    {
        bitpackPel = 0;

        //pixel 1 in bitpack image
        pel1=grayImage[yoff+bx];
        pel2=grayImage[yoff+bx+1];
        pel3=grayImage[yoff+bx+2];
        pel4=grayImage[yoff+bx+3];
        pel5=grayImage[yoff+bx+4];
        pel6=grayImage[yoff+bx+5];
        pel7=grayImage[yoff+bx+6];
        pel8=grayImage[yoff+bx+7];

        bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx]);
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+1] );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+2] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+3] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+4] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+5] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+6] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+7] );

        bitPackImage[byoff+(xOffsetLUT[bx])] = bitpackPel;

        //pixel 2 in bitpack image
        pel1=grayImage[yoff+bx+8];
        pel2=grayImage[yoff+bx+9];
        pel3=grayImage[yoff+bx+10];
        pel4=grayImage[yoff+bx+11];
        pel5=grayImage[yoff+bx+12];
        pel6=grayImage[yoff+bx+13];
        pel7=grayImage[yoff+bx+14];
        pel8=grayImage[yoff+bx+15];

        bitpackPel = 0;

        bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+8]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+9]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+10] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+11] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+12] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+13] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+14] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+15] );

        bitPackImage[byoff+(xOffsetLUT[bx+8])] = bitpackPel;

        //pixel 3 in bitpack image
        pel1=grayImage[yoff+bx+16];
        pel2=grayImage[yoff+bx+17];
        pel3=grayImage[yoff+bx+18];
        pel4=grayImage[yoff+bx+19];
        pel5=grayImage[yoff+bx+20];
        pel6=grayImage[yoff+bx+21];
        pel7=grayImage[yoff+bx+22];
        pel8=grayImage[yoff+bx+23];

        bitpackPel = 0;

        bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+16]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+17]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+18] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+19] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+20] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+21] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+22] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+23] );

        bitPackImage[byoff+(xOffsetLUT[bx+16])] = bitpackPel;

        //pixel 4 in bitpack image
        pel1=grayImage[yoff+bx+24];
        pel2=grayImage[yoff+bx+25];
        pel3=grayImage[yoff+bx+26];
        pel4=grayImage[yoff+bx+27];
        pel5=grayImage[yoff+bx+28];
        pel6=grayImage[yoff+bx+29];
        pel7=grayImage[yoff+bx+30];
        pel8=grayImage[yoff+bx+31];

        bitpackPel = 0;

        bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+24]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+25]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+26] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+27] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+28] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+29] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+30] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+31] );

        bitPackImage[byoff+(xOffsetLUT[bx+24])] = bitpackPel;
    }
}
}

I know this algorithm will potentially miss some trailing pixels in each row, but don't worry about that.

As you can see for every monochrome byte, I process 8 grayscale pixels.

Where you see pel8<=threshold is a neat little trick that resolves to 0 or 1 and is much faster than if{} else{}

For every increment of X I pack a bit into a higher order bit than the previous X

so for the first set of 8 pixels in the grayscale image

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8

This is what the bits in the byte look like (obviously each numbered bit is just the threshold result of having processed the corresponding numbered pixel but you get the idea)

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8

PHEW that should be it. Feel free to have some fun with some nifty bit twiddling tricks that would squeeze more juice out of this algorithm.

With compiler optimizations turned on this function takes on average 16 milliseconds on a roughly 5000 x 2200 pixel image on a core 2 duo machine.

EDIT:

R..'s suggestion was to remove the shift LUT and just use constants which is actually perfectly logical...I have modified the OR'ing of each pixel to be as such:

void binarize( unsigned char grayImage[], unsigned char bitPackImage[], int threshold, int grayscaleWidth, int grayscaleHeight, int  bitPackedScanlineStride, int bitpackedLength,  int shiftLUT[], int xOffsetLUT[] )
{
int yoff;
int byoff;
unsigned char bitpackPel=0;
unsigned char pel1=0;
unsigned char  pel2=0;
unsigned char  pel3=0;
unsigned char  pel4=0;
unsigned char  pel5=0;
unsigned char  pel6=0;
unsigned char  pel7=0;
unsigned char  pel8=0;
int checkX=grayscaleWidth-32;
int checkY=grayscaleHeight;

for ( int by = 0 ; by < checkY; by++)
{
    yoff=by*grayscaleWidth;
    byoff=by*bitPackedScanlineStride;

    for( int bx = 0; bx < checkX; bx+=32)
    {
        bitpackPel = 0;

        //pixel 1 in bitpack image
        pel1=grayImage[yoff+bx];
        pel2=grayImage[yoff+bx+1];
        pel3=grayImage[yoff+bx+2];
        pel4=grayImage[yoff+bx+3];
        pel5=grayImage[yoff+bx+4];
        pel6=grayImage[yoff+bx+5];
        pel7=grayImage[yoff+bx+6];
        pel8=grayImage[yoff+bx+7];

        /*bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx]);
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+1] );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+2] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+3] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+4] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+5] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+6] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+7] );*/
        bitpackPel |= ( (pel1<=threshold) << 7);
        bitpackPel |= ( (pel2<=threshold) << 6 );
        bitpackPel |= ( (pel3<=threshold) << 5 );
        bitpackPel |= ( (pel4<=threshold) << 4 );
        bitpackPel |= ( (pel5<=threshold) << 3 );
        bitpackPel |= ( (pel6<=threshold) << 2 );
        bitpackPel |= ( (pel7<=threshold) << 1 );
        bitpackPel |= ( (pel8<=threshold)  );

        bitPackImage[byoff+(xOffsetLUT[bx])] = bitpackPel;

        //pixel 2 in bitpack image
        pel1=grayImage[yoff+bx+8];
        pel2=grayImage[yoff+bx+9];
        pel3=grayImage[yoff+bx+10];
        pel4=grayImage[yoff+bx+11];
        pel5=grayImage[yoff+bx+12];
        pel6=grayImage[yoff+bx+13];
        pel7=grayImage[yoff+bx+14];
        pel8=grayImage[yoff+bx+15];

        bitpackPel = 0;

        /*bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+8]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+9]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+10] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+11] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+12] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+13] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+14] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+15] );*/
         bitpackPel |= ( (pel1<=threshold) << 7);
        bitpackPel |= ( (pel2<=threshold) << 6 );
        bitpackPel |= ( (pel3<=threshold) << 5 );
        bitpackPel |= ( (pel4<=threshold) << 4 );
        bitpackPel |= ( (pel5<=threshold) << 3 );
        bitpackPel |= ( (pel6<=threshold) << 2 );
        bitpackPel |= ( (pel7<=threshold) << 1 );
        bitpackPel |= ( (pel8<=threshold)  );


        bitPackImage[byoff+(xOffsetLUT[bx+8])] = bitpackPel;

        //pixel 3 in bitpack image
        pel1=grayImage[yoff+bx+16];
        pel2=grayImage[yoff+bx+17];
        pel3=grayImage[yoff+bx+18];
        pel4=grayImage[yoff+bx+19];
        pel5=grayImage[yoff+bx+20];
        pel6=grayImage[yoff+bx+21];
        pel7=grayImage[yoff+bx+22];
        pel8=grayImage[yoff+bx+23];

        bitpackPel = 0;

        /*bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+16]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+17]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+18] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+19] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+20] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+21] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+22] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+23] );*/
          bitpackPel |= ( (pel1<=threshold) << 7);
        bitpackPel |= ( (pel2<=threshold) << 6 );
        bitpackPel |= ( (pel3<=threshold) << 5 );
        bitpackPel |= ( (pel4<=threshold) << 4 );
        bitpackPel |= ( (pel5<=threshold) << 3 );
        bitpackPel |= ( (pel6<=threshold) << 2 );
        bitpackPel |= ( (pel7<=threshold) << 1 );
        bitpackPel |= ( (pel8<=threshold)  );


        bitPackImage[byoff+(xOffsetLUT[bx+16])] = bitpackPel;

        //pixel 4 in bitpack image
        pel1=grayImage[yoff+bx+24];
        pel2=grayImage[yoff+bx+25];
        pel3=grayImage[yoff+bx+26];
        pel4=grayImage[yoff+bx+27];
        pel5=grayImage[yoff+bx+28];
        pel6=grayImage[yoff+bx+29];
        pel7=grayImage[yoff+bx+30];
        pel8=grayImage[yoff+bx+31];

        bitpackPel = 0;

        /*bitpackPel |= ( (pel1<=threshold) << shiftLUT[bx+24]  );
        bitpackPel |= ( (pel2<=threshold) << shiftLUT[bx+25]  );
        bitpackPel |= ( (pel3<=threshold) << shiftLUT[bx+26] );
        bitpackPel |= ( (pel4<=threshold) << shiftLUT[bx+27] );
        bitpackPel |= ( (pel5<=threshold) << shiftLUT[bx+28] );
        bitpackPel |= ( (pel6<=threshold) << shiftLUT[bx+29] );
        bitpackPel |= ( (pel7<=threshold) << shiftLUT[bx+30] );
        bitpackPel |= ( (pel8<=threshold) << shiftLUT[bx+31] );*/
         bitpackPel |= ( (pel1<=threshold) << 7);
        bitpackPel |= ( (pel2<=threshold) << 6 );
        bitpackPel |= ( (pel3<=threshold) << 5 );
        bitpackPel |= ( (pel4<=threshold) << 4 );
        bitpackPel |= ( (pel5<=threshold) << 3 );
        bitpackPel |= ( (pel6<=threshold) << 2 );
        bitpackPel |= ( (pel7<=threshold) << 1 );
        bitpackPel |= ( (pel8<=threshold)  );


        bitPackImage[byoff+(xOffsetLUT[bx+24])] = bitpackPel;
    }
}
}

I am now testing on an Intel Xeon 5670, using (GCC) 4.1.2. Under these specs the hardcoded bitshift are 4 ms slower than using my original LUT algorithm. In the Xeon and GCC, the LUT algorithm takes on average 8.61 ms and the hard coded bitshift takes on average 12.285 ms.

Your lookup table is useless. Simply computing the shift (if you do it correctly, and not using the `%` operator with a signed integer, which is wastefully slow) is a lot faster than the table lookup. Or, much better, you could just unroll the loop and hard-code the 8 shifts. Typically bitshift by a constant is a lot faster than variable bitshift, so it could help a lot. — R.. GitHub STOP HELPING ICE, Sep 14 '10 at 00:39
I have modified the algorithm to simply use constant bitshifts.... it actually ends up being 4 ms slower than the LUT. I'm now using an Intel Xeon on GCC 1.4.2. The algorithm using the LUT takes 8.61 ms on average and without the LUT it takes on average 12.285 ms. — alessandro ferrucci, Sep 14 '10 at 01:25
+1 to R.., and the second lut is equally useless as `x/8` will become `x>>3`, which is faster than `*(lut+x)`, since you dont need to dereference a pointer. If you really think that insane portability is worthwhile (and not precluded by other constructs you're using) then you can use `x/CHAR_BIT`. — David X, Sep 14 '10 at 01:29
@alssandro, that doesnt sound right, can you post the code you used to get 8.61 vs 12.285? — David X, Sep 14 '10 at 01:32
`x/8` will not become `x>>3` unless `x` is unsigned or the compiler can determine that `x` will never be negative. — R.. GitHub STOP HELPING ICE, Sep 14 '10 at 03:22
BTW, this entire function should be a few lines long, not the monstrosity you've made it... See the approach I gave in an answer. — R.. GitHub STOP HELPING ICE, Sep 14 '10 at 03:24
@David X: Note I have edited the question to include the function where the shift does not include a LUT. — alessandro ferrucci, Sep 14 '10 at 10:18
@R, You are right the function shouldn't be that long, the only reason it was this long is that I just tried unrolling the loop to process 32 packed pixels at a time. — alessandro ferrucci, Sep 14 '10 at 10:21
@alessandro, okay, I *think* the speedup is because GCC is using SIMD instructions to do the `<=` and/or `<<`, but it should be able to do that with constant values too. Could you post the assemby (`gcc -S`)? (Note: if you dont want to make your question longer, you could try a code hosting site like `codepad.org`, `ideone.com`, or `pastebin.com`.) — David X, Sep 14 '10 at 23:51

score 2 · Accepted Answer · answered Sep 14 '10 at 01:49

2

Try something like:

unsigned i, w8=w>>3, x;
for (i=0; i<w8; i++) {
    x = thres-src[0]>>1&0x80;
    x |= thres-src[1]>>2&0x40;
    x |= thres-src[2]>>3&0x20;
    x |= thres-src[3]>>4&0x10;
    x |= thres-src[4]>>5&0x08;
    x |= thres-src[5]>>6&0x04;
    x |= thres-src[6]>>7&0x02;
    x |= thres-src[7]>>8&0x01;
    out[i] = x;
    src += 8;
}

You can figure out the extra code for the remainder at the end of the line of the width is not a multiple of 8, or you could just pad/align the source to ensure it's a multiple of 8.

answered Sep 14 '10 at 01:49

R.. GitHub STOP HELPING ICE

208,859
35
376
711

1

Are you sure those shifts couldn't go from 0 to 7, rather than 1 to 8 (assuming that the threshold and src are both 8 bit values). – caf Sep 14 '10 at 04:50
Yes, I chose the correct shift values. I'm shifting down bit 8, not bit 7, because I want the borrow bit from the integer result. Bit 7 could be 0 or 1 regardless of whether `thres-src[k]` wrapped modulo `UINT_MAX+1`. – R.. GitHub STOP HELPING ICE Sep 14 '10 at 12:36

Paul R · Answer 2 · 2013-07-19T06:11:34.657

You can do this with SSE quite easily, processing 16 pixels at a time, e.g.

load vector (16 x 8 bit unsigned)
add (255 - threshold) to each element
use PMOVMSKB to extract sign bits into 16 bit word
store 16 bit word

Example code using SSE intrinsics (warning: untested !):

void threshold_and_pack(
    const uint8_t * in_image,       // input image, 16 byte aligned, height rows x width cols, width = multiple of 16
    uint8_t * out_image,            // output image, 2 byte aligned, height rows x width/8 cols, width = multiple of 2
    const uint8_t threshold,        // threshold
    const int width,
    const int height)
{
    const __m128i vThreshold = _mm_set1_epi8(255 - threshold);
    int i, j;

    for (i = 0; i < height; ++i)
    {
        const __m128i * p_in = (__m128i *)&in_image[i * width];
        uint16_t * p_out = (uint16_t *)&out_image[i * width / CHAR_BIT];

        for (j = 0; j < width; j += 16)
        {
            __m128i v = _mm_load_si128(p_in);
            uint16_t b;

            v = _mm_add_epi8(v, vThreshold);
            b = _mm_movemask_epi8(v);   // use PMOVMSKB to pack sign bits into 16 bit word

            *p_out = b;

            p_in++;
            p_out++;
        }
    }
}

fast threshold and bit packing algorithm ( possible improvements ? )

2 Answers2

Linked