30 September 2012

Python code for non rigid registration using opencv and numpy

Here is a simple python code for registering two images. The transform is non-rigid. The error metric is sum of squared difference. Code is very simple, does not handle physical spacing. Main aim was to get hands dirty of numpy/scipy, python, opencv cv2 interface etc.
One practical tip for registration: if all else fails, smooth heavily the motion fields ;-)

##################################################
"""
Created on Sun Jan 15 10:41:19 2012

@author: vsk
"""
import sys, os, cv2, numpy, matplotlib, pylab
def SSD_Register(static_image, moving_image):
    IMG_BLUR = 11;
    ITER = 200;
    TIME_STEP = 0.00010;
    SMOOTH = 1.0;
    K_SIZE = 11;
   
    # Create float versions of all images
    registered_image = numpy.empty(static_image.shape, dtype=numpy.float32);  
    static_image_float = numpy.empty(static_image.shape, dtype=numpy.float32);
    mov_image_float = numpy.empty(static_image.shape, dtype=numpy.float32);

    for i in range(static_image.shape[0]):
        for j in range(static_image.shape[1]):
            static_image_float[i,j] = static_image[i,j];
            mov_image_float[i,j] = moving_image[i,j];
                       
    static_image_blur = cv2.GaussianBlur(static_image_float,(IMG_BLUR,IMG_BLUR), 1.20);
    mov_image_float = cv2.GaussianBlur(mov_image_float, (IMG_BLUR,IMG_BLUR), 1.20);
   
    difference_image = numpy.empty(static_image.shape, dtype = numpy.float32);
    # Space for motion field
    x_field = numpy.empty(static_image.shape, dtype=numpy.float32);
    y_field = numpy.empty(static_image.shape, dtype=numpy.float32);
    for i in range(x_field.shape[0]):
        for j in range(x_field.shape[1]):
            x_field[i,j] = j;
            y_field[i,j] = i;
           
    x_field_update = numpy.empty_like(x_field);
    y_field_update = numpy.empty_like(y_field);
    # Gradient image
    grad_x = numpy.empty_like(x_field_update);
    grad_y = numpy.empty_like(x_field_update);
   
    grad_x = cv2.Sobel(mov_image_float, cv2.CV_32FC1, 1, 0);
    grad_y = cv2.Sobel(mov_image_float, cv2.CV_32FC1, 0, 1);
    grad_x_resampled = numpy.copy(grad_x);
    grad_y_resampled = numpy.copy(grad_y);
    # registration iterations begin
   
    for i in range(x_field.shape[0]):
        for j in range(x_field.shape[1]):
            registered_image[i,j] = mov_image_float[i,j];
           
    difference_image = registered_image - static_image_float;
    print(numpy.sum(numpy.abs(difference_image))),
    ssd_err = numpy.zeros((ITER,1));
    for i in range(ITER):
    #    print i,
        difference_image = registered_image - static_image_float;
        #print(numpy.max(difference_image)),
        #print(numpy.min(difference_image))
    #    print(numpy.sum(numpy.abs(difference_image))),
        ssd_err[i]= numpy.sum(numpy.abs(difference_image));
        x_field_update = TIME_STEP*difference_image*grad_x_resampled;
        y_field_update = TIME_STEP*difference_image*grad_y_resampled;

        x_field_update = cv2.GaussianBlur(x_field_update, (K_SIZE, K_SIZE), 2.5);
        y_field_update = cv2.GaussianBlur(y_field_update, (K_SIZE, K_SIZE), 2.5);
    #    print(numpy.max(y_field_update))
        x_field = x_field - x_field_update;
        y_field = y_field - y_field_update;

        registered_image = cv2.remap(mov_image_float, x_field, y_field, cv2.INTER_CUBIC);
        grad_x_resampled = cv2.remap(grad_x, x_field, y_field, cv2.INTER_CUBIC);
        grad_y_resampled = cv2.remap(grad_y, x_field, y_field, cv2.INTER_CUBIC);
    print('')
    #pylab.imshow(registered_image); pylab.show(); pylab.gray()
    print(numpy.sum(numpy.abs(difference_image)))
    # display section
    display_image = numpy.copy(static_image);
    for i in range(x_field.shape[0]):
        for j in range(x_field.shape[1]):
            display_image[i,j] = registered_image[i,j];

#     cv2.namedWindow('Registered');
#     cv2.imshow("Registered", display_image);
#    
#     cv2.namedWindow('Final');
#     cv2.imshow("Final", moving_image);
#
#     cv2.namedWindow('Initial');
#     cv2.imshow("Initial", static_image);
#    
#     cv2.waitKey(-400)
#     cv2.destroyAllWindows();
#     pylab.plot(ssd_err)
#     pylab.show()
   
    pylab.subplot(131)
    pylab.imshow(static_image)
    frame = pylab.gca();
    frame.axes.get_xaxis().set_visible( False);
    frame.axes.get_yaxis().set_visible( False);
    pylab.title('Static Image');
   
    pylab.subplot(133)
    pylab.imshow(moving_image)
    frame = pylab.gca();
    frame.axes.get_xaxis().set_visible( False);
    frame.axes.get_yaxis().set_visible( False);
    pylab.title('Moving Image');
   
    pylab.subplot(132)
    pylab.imshow(registered_image),
    frame = pylab.gca();
    frame.axes.get_xaxis().set_visible( False);
    frame.axes.get_yaxis().set_visible( False);
    pylab.title('Registered Image');
   
    pylab.show()
    return (registered_image, ssd_err);
   
# ================  Read system arguments else generate images =================
if len(sys.argv) != 3:
    #print('Usage: static_image moving_image');
    static_im = numpy.zeros( (256, 256), dtype=numpy.ubyte);
    static_im[100:150, 100:150] = 100;

    x_field = numpy.empty(static_im.shape, dtype=numpy.float32);
    y_field = numpy.empty(static_im.shape, dtype=numpy.float32);
    for i in range(x_field.shape[0]):
        for j in range(x_field.shape[1]):
            x_field[i,j] = j + 12.00485;
            y_field[i,j] = i -  8.35;
   
    mov_im = numpy.zeros( (256, 256), dtype=numpy.ubyte);
    mov_im = cv2.remap(static_im, x_field, y_field, cv2.INTER_CUBIC);
    #mov_im[102:152, 100:150] = 100;
   
    SSD_Register(static_im, mov_im);  
else:
    static_image = cv2.imread(sys.argv[1], cv2.CV_LOAD_IMAGE_GRAYSCALE);
    moving_image = cv2.imread(sys.argv[2], cv2.CV_LOAD_IMAGE_GRAYSCALE)
    SSD_Register(static_image, moving_image)

10 February 2012

Regarding a trip to Puri temple

I have been meaning to write this since a long time, but never got around to do so. I visited the Jagannath Puri temple in 2008. It is 4 years now, but the visit remains extremely fresh in my mind. I had heard of crowds at Puri from relatives and friends before. The rudeness, greed and filthy language used by pandas was also indicated at. The advice given to me or rather us(since we were travelling in a group), was that avoid pandas at all costs, even if you have to listen to their entreaties which gradually become threats and end in abuse. This happens especially if there are females in the group, and fortunately or otherwise there was one in our group. As soon as we got down from the vehicle, we were surrounded. There was one guy who said who promised to take a sum which for his guidance. We made it clear that we are not interested in making some special puja(specially priced, if one may add). He agreed but obviously after entering premises, he said that prasad needs to be "bought", else visit is not possible. Needless to say, though we were not interested and wanted to get rid of him, the female of the group, who got a bit emotional with all the flowery language, agreed to buy some 150/- "basket". The puja coupon counter is another mess, manned by a fat scoundrel. Be as it may, after waiting for quite sometime, in which the scoundrel rudely asked the lady to pipe down, the panda came back with a basket. After entering the main temple, he asked each of us to put money as it was some custom. At this point we had enough and asked him to go away, which he did after muttering many things unpleasant. As I was moving along the queue for darshan, the panda thief inside, waved a large basket in front of me and shouted or demanded to put something inside. I said nothing but contented by just glaring at him. He kind of piped down. All in all, it was a filthy experience. 
It is sad that a temple and a holy place, where thousands come to pray, seek mental peace, succour, and relief, give their thanks, are treated so shabbily, abused, shoved and made to undergo such pain. Also, it is all the more surprising because the Puri temple's administrative board has some King as a member. Clearly, someone is snoozing. It is high time that the pandas be given the same treatment as thugees were given by the British. One may think that I am over-reacting, but I have decided not enter the temple again until I hear that things have improved. The pandas should realise that they are playing with fire in treating the devotees in this fashion. It is sad that no one is Odisha is interested in improving the situation.

Mahabhoj

मन्नू भंडारी द्वारा लिखित "महाभोज" पढ़ रहा हूँ। बहुत मोटी किताब नहीं है। कहानी का विषय है राजनीति। सारांश कुछ ऐसी है। किसी एक छोटे से गाँव में एक हरिजन समाजसेवक का लाश मिलने पर सरकार और विपक्षी नेता, दोनों इस "मौके" का फायदा उठाना चाहते हैं। कारण है कि इस घटना के कुछ ही दिनों बाद चुनाव है। वैसे तो है मामूली चुनाव लेकिन विपक्षी नेता इसी से अपना "नव निर्माण" करना चाहते हैं, और ये भी स्वाभाविक है की सरकारी दल, जिसमे कुछ नेता असंतुष्ट हैं,  इतनी आसानी से हारने वाली नहीं है। अब तक तो आधी ही पढ़ी है, लेकिन इतना तो समझ में आ ही गया है कि चुनाव में कोई भी जीते, असली में हारने वाले तो गाँव वाले ही हैं। इसलिए कह रहा हूँ कि कहानी बहुत ही सहज और स्वाभाविक रूप में लिखी गयी है। विभिन्न नेताओं का ऐसा सहज चरित्रांकन, बिना एक शब्द व्यर्थ किये, बहुत अच्छा लगा। ऐसा नहीं लगता है कि ये कोई कहानी है, बल्कि राजनीती के भीतरी एवं बाहरी चालबाज़ी का आँखों देखा हाल है।
खरीदिये और पढ़िए!

14 January 2012

Sony Ericsson Elm quick review

This is a quick review of the Sony Elm mobile cellphone.
Motivation and Decision!
The reason for the purchase was that I was looking for a cellphone within Rs. 8,000-10,000, with reasonable camera performance. I was tempted by LG and especially Samsung offerings. Since I have a very positive experience with Samsung E2232, I had decided upon Samsung Galaxy Y, which is a cheaper version of the more pricey ones of the same name. I was not specifically looking for Android phones or touchscreen phones. I am slightly wary of the touchscreen technology coming at too low a price; I guess that any company has to maintain its margin and will cut down on quality if pushed too hard. Now, all the reviews for Samsung were great except for a few which said that the camera was barely passable. I was also not looking for high megapixel, 3 would be more than enough. Sadly, moronic customers push the companies toward megapixel wars leading to degradation of images, or maybe the companies attract morons by such number games. However, I hesitated after reading the poor quality images. Then, I chanced upon the Sony Ericsson Elm. On  paper it seemed great; a Sony Ericsson product with decent reviews. I also saw some sample images at a website. Seemed reasonable to me for a phone costing 8000 or so. So I bought it.

General Comments
As far as looks go, this is fairly good. It is very conventional, no scratching or touching anywhere. The keypads are good, I played some games(pre-loaded) without worrying about damaging the pad. I would say that, this is about the only thing which I have positive to say about this overpriced brick. The phone itself is very slow. I never thought I would pay much attention to the themes, but the utter dullness and the available themes is galling. The sound quality is miserable, or to be honest maybe i expected better, since i was not paying for touchscreen and Android, I thought that I was getting better hardware.
Photos!
Ah, this is what I bought the phone for?! 

No, certainly not! This is bright sunshine! The hopeless quality is immediately apparent. Based on this short experience of 1 week, I think that I made a wrong choice. My general suggestion is that if you want a camera, buy one and if you want a phone, buy one with such features as you wish. Do not mix and match both, at least not for anything costing less than 25,000. At that price or a bit more, you may get an entry level DSLR nowadays. However, if you do feel like donating money by spending it on such a combo, while I do not know which to recommend, I most certainly which NOT to go for: Sony Ericsson Elm!