Tuesday, March 8, 2016

How to devise an auto-encoder

I think Autoencoder is an important type of MLP since one way it's to some extent an extension of the traditional classification method utilizing MLP, another way it's the key to understanding some ingenious training ideal which weights high in application of Deep Learning.

In this post I will guide to build a very simple autoencoder from implementation perspective. I have to state first that I am not intended to break Ng's law that do not share source code, typically when most content of this post is borrowed from his famous website http://ufldl.stanford.edu/wiki/index.php/UFLDL_Tutorial. This problem is most of the time, even I honestly write the code myself, but some time later I have no idea of my original thinking. So I write here mainly for reference. The encouragement is be familiar with it, not directly copy it.

The implementation is in Matlab script, typically, the work needs to be done in the sparseAutoencoderCost.m. I would think I have had implemented the autoencoder correctly, however no guarantee:

function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ...
                                             lambda, sparsityParam, beta, data)

% visibleSize: the number of input units (probably 64)
% hiddenSize: the number of hidden units (probably 25)
% lambda: weight decay parameter
% sparsityParam: The desired average activation for the hidden units (denoted in the lecture
%                           notes by the greek alphabet rho, which looks like a lower-case "p").
% beta: weight of sparsity penalty term
% data: Our 64x10000 matrix containing the training data.  So, data(:,i) is the i-th training example.
 
% The input theta is a vector (because minFunc expects the parameters to be a vector).
% We first convert theta to the (W1, W2, b1, b2) matrix/vector format, so that this
% follows the notation convention of the lecture notes.

W1 = reshape(theta(1:hiddenSize*visibleSize), hiddenSize, visibleSize);
W2 = reshape(theta(hiddenSize*visibleSize+1:2*hiddenSize*visibleSize), visibleSize, hiddenSize);
b1 = theta(2*hiddenSize*visibleSize+1:2*hiddenSize*visibleSize+hiddenSize);
b2 = theta(2*hiddenSize*visibleSize+hiddenSize+1:end);

% Cost and gradient variables (your code needs to compute these values).
% Here, we initialize them to zeros.
cost = 0;
W1grad = zeros(size(W1));
W2grad = zeros(size(W2));
b1grad = zeros(size(b1));
b2grad = zeros(size(b2));

%% ---------- YOUR CODE HERE --------------------------------------
%  Instructions: Compute the cost/optimization objective J_sparse(W,b) for the Sparse Autoencoder,
%                and the corresponding gradients W1grad, W2grad, b1grad, b2grad.
%
% W1grad, W2grad, b1grad and b2grad should be computed using backpropagation.
% Note that W1grad has the same dimensions as W1, b1grad has the same dimensions
% as b1, etc.  Your code should set W1grad to be the partial derivative of J_sparse(W,b) with
% respect to W1.  I.e., W1grad(i,j) should be the partial derivative of J_sparse(W,b)
% with respect to the input parameter W1(i,j).  Thus, W1grad should be equal to the term
% [(1/m) \Delta W^{(1)} + \lambda W^{(1)}] in the last block of pseudo-code in Section 2.2
% of the lecture notes (and similarly for W2grad, b1grad, b2grad).
%
% Stated differently, if we were using batch gradient descent to optimize the parameters,
% the gradient descent update to W1 would be W1 := W1 - alpha * W1grad, and similarly for W2, b1, b2.
%

% Since it's a special case, so iteration is avoided here
% first we calculate the activation of the first hidden layer

% Since dimension of W1 is of hiddenSize*visibleSize, so W'*x is a proposed
% alternative, since the bias term needs to be taken into consideration, we
% modify W to reflect such a effectiveness

% Keep in mind it's intended to utilize the computation in vector form as thorough as possible

ss = size(data, 2);

% Catering for the bias term, so we manipulate the weights and inputs 
% for the input layer
W1_e = [b1 W1];
data_e = [ones(1, ss); data];

% Calculation for the hidden layer
u1 = W1_e * data_e;
x2 = sigmoid(u1);


% The same manipulation applies to the input layer
W2_e = [b2 W2];
x2_e = [ones(1, ss); x2];


% Calculation for the output layer
u2 = W2_e * x2_e;
output = sigmoid(u2);

% Up to now, the forward propagation is completed

% Now it's time to build the cost function

% First the avarage cost without weight decay could be calculated
diff_output = output - data;

cost = mean(sum(diff_output .^ 2)) / 2;

% Take the weight decay into consideration
cost = cost + lambda * (sum(sum(W1 .^ 2)) + sum(sum(W2 .^ 2))) / 2;

% The average activity of the hidden layer can be computed as follows
rho = mean(x2, 2);

sparsity_cost = sparsityParam * log(sparsityParam ./ rho) + ...
    (1 - sparsityParam) * log((1 - sparsityParam) ./ (1 - rho));

cost = cost + beta * sum(sparsity_cost);

% Now we compute the gradient via back-propagation
% We do the computation from an element-wise viewpoint but from an
% vector-wise way
delta_output = (output - data) .* output .* (1 - output);

b2grad = mean(delta_output, 2);

for i = 1 : ss;
    W2grad = W2grad + delta_output(:, i) * x2(:, i)';
end

W2grad = W2grad / ss + lambda * W2;

delta_sparsity = -1 * sparsityParam ./ rho + ...
    (1 - sparsityParam) ./ (1 - rho);

delta_hidden = W2' * delta_output + repmat(beta * delta_sparsity, 1, ss);

delta_hidden = delta_hidden .* x2 .* (1 - x2);

b1grad = mean (delta_hidden, 2);

for i = 1 : ss;
    W1grad = W1grad + delta_hidden(:, i) * data(:, i)';
end

W1grad = W1grad / ss + lambda * W1;


%-------------------------------------------------------------------
% After computing the cost and gradient, we will convert the gradients back
% to a vector format (suitable for minFunc).  Specifically, we will unroll
% your gradient matrices into a vector.

grad = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)];

end

%-------------------------------------------------------------------
% Here's an implementation of the sigmoid function, which you may find useful
% in your computation of the costs and the gradients.  This inputs a (row or
% column) vector (say (z1, z2, z3)) and returns (f(z1), f(z2), f(z3)). 

function sigm = sigmoid(x)
  
    sigm = 1 ./ (1 + exp(-x));
end


It's estimated the most difficulty part is the implementation of Autoencoder, once we have one, some interesting results can be spotted. The first figure belows show some samples to be learnt from, and the second figure is the weights learnt however displayed in an squared image way:



Happy exploring, happy coding.

Friday, December 25, 2015

How to draw OpenCV Mat data into native Windows client area utilizing GDI


under some scenario, probably one wants to combine the OpenCV's ability to manipulate pixel values, and the Windows ability to display the resulting image in a native way (for example, utilizing the GDI system). Such a requirement demands some deep understanding of image representation between these two systems, and probably a non-trivial task than just taking for granted.


I came across a same situation as just mentioned above, and after some try-and-error iteration, finally everything seems have been put onto track. since it hasn't undergone thorough test, I am not sure the tracks are all right tracks.

would like to share here for polishing or beneficial for others facing the same situation.

The header file is as follows:

#pragma once

#include <opencv2/core.hpp>

class COpenCVImage
{
public:
COpenCVImage(LPCTSTR lpszFileName);
~COpenCVImage(void);

BOOL GetSize(LPRECT lpRect);
BOOL DrawBitmap(HDC hdc, LPRECT lpRect);

BOOL Sharpen(HWND hWnd);
BOOL Flip(HWND hWnd);

private:
bool createBitmap(HDC hdc, cv::Mat mat);
private:
cv::Mat m_matImage;
HBITMAP m_hBitmap;

bool m_bInit;
};


the implementation file is as follows:

#include "stdafx.h"
#include "OpenCVImage.h"

#include <Shlwapi.h>
// The following header file will be excluded, since we will draw the picture 
// to the window all by ourselves
// #include <opencv2/highgui.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>

COpenCVImage::COpenCVImage(LPCTSTR lpszFileName) : m_hBitmap(NULL), m_bInit(false)
{
if(!PathFileExists(lpszFileName))
{
MessageBox(NULL, _T("Cannot open the image file"), _T("Error"), MB_OK);
return ;
}

char* pFileName;
UINT nLen;

nLen = WideCharToMultiByte(CP_UTF8, 0, lpszFileName, sizeof(TCHAR) * (_tcslen(lpszFileName) + 1), NULL, 0, NULL, NULL);
pFileName = new char[nLen];
WideCharToMultiByte(CP_UTF8, 0, lpszFileName, sizeof(TCHAR) * (_tcslen(lpszFileName) + 1), pFileName, nLen, NULL, NULL);

m_matImage = cv::imread(pFileName, cv::IMREAD_COLOR);
delete [] pFileName;

if (m_matImage.empty())
{
MessageBox(NULL, _T("Cannot read in the image file"), _T("Error"), MB_OK);
return;
}

m_bInit = true;
}

COpenCVImage::~COpenCVImage(void)
{
if (m_hBitmap)
DeleteObject(m_hBitmap);
}


BOOL COpenCVImage::DrawBitmap(HDC hdc, LPRECT lpRect)
{
if (!m_bInit)
return FALSE;

HDC hDCMem = CreateCompatibleDC(hdc);

if (!m_hBitmap)
{
createBitmap(hDCMem, m_matImage);
}

SelectObject(hDCMem, m_hBitmap);

BitBlt(hdc, lpRect->left, lpRect->top, lpRect->right - lpRect->left, 
lpRect->bottom - lpRect->top, hDCMem, 0, 0, SRCCOPY);

DeleteDC(hDCMem);

return TRUE;
}


BOOL COpenCVImage::GetSize(LPRECT lpRect)
{
if (!m_bInit)
{
memset(lpRect, 0, sizeof(RECT));
return FALSE;
}

lpRect->left = 0;
lpRect->top = 0;
lpRect->right = m_matImage.cols;
lpRect->bottom = m_matImage.rows;

return TRUE;
}

BOOL COpenCVImage::Sharpen(HWND hWnd)
{
if (!m_bInit)
return FALSE;

cv::Mat matFiltImage;
cv::Matx33f matKernel;

matKernel(1, 1) = 5.0;
matKernel(0, 1) = - 1.0;
matKernel(1, 0) = - 1.0;
matKernel(1, 2) = - 1.0;
matKernel(2, 1) = - 1.0;

cv::filter2D(m_matImage, matFiltImage, m_matImage.depth(), matKernel);

HDC hDCMem = CreateCompatibleDC(GetDC(hWnd));
createBitmap(hDCMem, matFiltImage);
DeleteObject(hDCMem);

InvalidateRect(hWnd, NULL, TRUE);

return TRUE;
}

BOOL COpenCVImage::Flip(HWND hWnd)
{
if (!m_bInit)
return FALSE;

cv::Mat matFlipImage;
cv::flip(m_matImage, matFlipImage, 1);

HDC hDCMem = CreateCompatibleDC(GetDC(hWnd));
createBitmap(hDCMem, matFlipImage);
DeleteObject(hDCMem);

InvalidateRect(hWnd, NULL, TRUE);

return TRUE;
}

bool COpenCVImage::createBitmap(HDC hdc, cv::Mat mat)
{
unsigned char* lpBitmapBits;

if (m_hBitmap)
DeleteObject(m_hBitmap);

BITMAPINFO bi; 
ZeroMemory(&bi, sizeof(BITMAPINFO));
bi.bmiHeader.biSize = sizeof(BITMAPINFOHEADER);
bi.bmiHeader.biWidth = m_matImage.cols;
bi.bmiHeader.biHeight = - m_matImage.rows;
bi.bmiHeader.biPlanes = 1;
bi.bmiHeader.biBitCount = (m_matImage.elemSize() << 3);
bi.bmiHeader.biCompression = BI_RGB;

m_hBitmap = CreateDIBSection(hdc, &bi, DIB_RGB_COLORS, (VOID**)&lpBitmapBits, NULL, 0);

#define ALIGN(x,a)              __ALIGN_MASK(x, a-1)
#define __ALIGN_MASK(x,mask)    (((x)+(mask))&~(mask))
int width = m_matImage.cols * m_matImage.elemSize();
int pitch = ALIGN(width, 4);
#undef __ALIGN_MASK
#undef ALIGN

for(int i = 0; i < m_matImage.rows; i ++)
{
unsigned char* data = mat.ptr<unsigned char>(i);
memcpy(lpBitmapBits, data, width);
lpBitmapBits += pitch;
}

return true;
}

a typical usage probably be:

1. instantiate the instance, for example:
COpenCVImage image(_T("lenna.jpg"));

2. respond to the WM_PAINT message:

case WM_PAINT:
hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code here...

RECT rect;
image.GetSize(&rect);
image.DrawBitmap(hdc, &rect);

EndPaint(hWnd, &ps);

break;


the final result is as follows:


appreciate any suggestion to enhance the code

Sunday, December 20, 2015

Fourier in image processing (I)

before indulge into the frequency field of image processing, two points must be manifested:

1. the mathematical rigor is a good thing, however in practice just stick to it in some extent in case of seeking trouble for oneself. so don't be panic if you see someone exchange summation and integration freely.

2. casualness doesn't give one enough reason to neglect some basic concepts concerning with Fourier theory, like the difference between Fourier series and Fourier transform. otherwise, your further deduction probably in the wrong track.

fft: fast Fourier transform for 1-D vector
Note, to really utilize the virtue of fft, the size of vector must a power of 2. otherwise it will fall back to the general discrete Fourier transform
example:
u = [1, 2, 4, 1];
F = fft(u);

the result is as follows:
F = [8.0000 + 0.0000i  -3.0000 - 1.0000i   2.0000 + 0.0000i  -3.0000 + 1.0000i]

fft2: fast Fourier transform for 2-D matrix
identical to fft, if the size of matrix along some dimension is not the power of 2, along that dimension, it automatically falls back into discrete Fourier transform.
example:
S = fspecial('sobel');
F = fft2(S);

the results are as follows:
S = [

     1     2     1
     0     0     0
    -1    -2    -1
];
F = [
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   6.0000 - 3.4641i  -1.5000 - 0.8660i  -0.0000 + 1.7321i
   6.0000 + 3.4641i  -0.0000 - 1.7321i  -1.5000 + 0.8660i
];

fftshift: shift the DC component to the center in order to benefit inspection

example:
FC = fftshift(F);

the result is as follows:
FC = [

  -1.5000 + 0.8660i   6.0000 + 3.4641i  -0.0000 - 1.7321i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  -0.0000 + 1.7321i   6.0000 - 3.4641i  -1.5000 - 0.8660i
];

ifftshift: inverse of fftshift

following, i will give a comprehensive example to demonstrate the usage of above functions. the purpose of the example is apply low-pass filter in the frequency domain to some image, so as to observe the results:

function freqfilt(img_path)

try
    img = imread(img_path);
catch
    msgbox('cannot open image file, check the file name!');
    return ;
end

h1 = figure(1); imshow(img);
figure(h1); title('Original image');

% Convert into gray-level image if necessary. you've been warned

% FFT transform
img_freq = fft2(img);

% Shift the DC component into center for observing
img_freq = fftshift(img_freq);

% Adjust the frequency value for displaying
img_freq_mag = abs(img_freq);
img_freq_mod = log(img_freq_mag + 1);

% Normalize the matrix since elements are of floating point datatype

f_max = max(max(img_freq_mod));
img_freq_mod = img_freq_mod / f_max;

h2 = figure(2); imshow(img_freq_mod);
figure(h2); title('Spectrum distribution');

% Create a mask for frequency filed filting
[cx, cy] = size(img);

mask = bsxfun(@plus, ((1:cx) - ceil(cx/2)) .^ 2, (transpose(1:cy) - ceil(cy/2)) .^ 2) < 100^2;
h3 = figure(3); imshow(uint8(mask * 255));
figure(h3); title('Mask image');

% Do the filting
img_freq = img_freq .* mask;

% Shift the freqency back to do inverse fft
img_freq = ifftshift(img_freq);

% Note the relation between ifft by utilizing forwarding fft
img_freq = img_freq';
img_ift = fft2(img_freq);
img_ift = img_ift' / numel(img);

% Trucate to real numbers to surpress round-off error
img_ift = real(img_ift);

h4 = figure(4); imshow(uint8(img_ift));
figure(h4); title('Filtered image');

the results are as follows:



imfilter: filtering a given image with specific kernel
example:
img = imfilter(rgb2gray(imread('test.jpg')), fspecial('guassian'));

note, for all kinds of filters, a normalization is recommended before applying to the image. for low pass filter, it's not necessary for adding the original image for displaying. generally for high pass filter, it's tended to request adding the original image before displaying.

freqz2: frequency response for 2-D matrices.
note, this function is mainly for testing kernels, architect your own utility to display the spectrum of some large image.
example:
freqz2(fspecial('laplacian'));

the result is as follows:

Friday, December 18, 2015

Matlab functions and image processing toolbox functions summary (I)

imread: load an image into matlab workspace

imshow: show an loaded image or synthetic image by some special matrix

example:
img_color = imread('forest.jpeg', 'jpg');
imshow(img_color);

the figure being shown is as follows:
rgb2gray: convert a color image to an gray-level image

example: 
img_gray = rgb2gray(img_color);


imhist: demonstrate the histogram of an image. note this routine can only apply to gray-level image. the answer is color image is generally abundant with intensity values, (255x255x255), such a sparse intensity over weight demonstration is of little meaningfulness.

example:
imhist(img_gray);

the result is as follows:
histeq: do the histogram equalization to specific image

example:
img_eq = histeq(img_gray);
imshow(img_eq);

the result is as follows:
graythresh: compute the threshold for a given image
im2bw: convert an image into black-white image

note above the above function are tended to use jointly
gt = graythresh(img_gray);
img_bw = im2bw(img_gray, gt);
imshow(img_bw);

the result is as follows:
note the above example is operated upon image with data type unsigned char. however it could also be operated on float number. refer to the following example, but before that, let's introduce another function:
im2double: convert the intensity value into double float number

so the final example is:
img_f = im2double(img_gray);
imshow(im2bw(img_f, graythresh(img_f)));

it will display the identical image.

imadjust: adjust the intensity level of an given image
the intensity level transformation tends to be governed by several common functions, like exponential functions, and piecewise linear functions. however since imadjust saturate the intensity value at two input scale ends, it's probably in helping inspect the image instead of operating on the image:

stretchlim: obtain the upper and lower bound of intensity value

example:
img_adj = imadjust(img_f, stretchlim(img_f), [0, 1], 0.25);
imshow(img_adj);
an explanation of the example. it will map the intensity value between the minimum and maximum intensity of img_f into the full range [0, 1]. here 0.25 means the exponential value of the transform function s=r^alpha, r is the old intensity value and s is the new intensity value.

the result is as follows:




sometimes, it's necessary obtain the negative image of the original one, it can be completed by the following way:
img_neg = imadjust(img_f, [0, 1], [1, 0]);

there's also a function fulfil the same purpose:
imcomplement: obtain the complement image of the original one
example:
img_neg = imcomplement(img_f);

the result is as follows:



Reference:
http://www.cs.otago.ac.nz/cosc451/Resources/matlab_ipt_tutorial.pdf
Digital Image Processing, 3e, Gonzalez & Woods






Special functions related to image processing

The following just summaries some common functions in Matlab related to image processing. it will be updated if new ones are chanced upon and worthy-while.

conv: convolution of two 1-dim vectors

example 1:
suppose u = [1, 2, 3]; v = [1, -1], then r = conv(u, v) = [1, 1, 1, -3]
the procedure is as follows:
1. flip v, so we have v = [-1, 1]
2. overlap the last element of v, that's 1, with the first element of u, that's 1, do the multiplication and addition. note if we assign u[0] = 1, u[1] = 2 and u[2] = 3, we implicitly padding u[-1] = 0 ahead of u. so we have u[-1] * 1 + u[0]*-1 = 0*(-1) + 1*1 = 1
3. advance v or slip v along the positive direction for one step, and do the multiplication and addition. so we have 1*(-1) + 2*1 = 1
4. repeat the above procedures, we will have the final result
note: conv(u, v) is equivalent to conv(u, v, 'full')

example 2:
suppose u = [1, 2, 3]; v = [1, -1], then r = conv(u, v, 'same') = [1, 2, -3];
the procedures is as follows:
1. flip v, so we have v = [-1, 1]
2. align u, v to the same head elements. alternatively, if we denote u[0] = 1, u[1] = 2, u[2] = 3, then we have v[0] = -1, v[1] = 1
3. do the multiplication and addition until the last chance of having u and v overlapping. that's means the final value must be u[end]*v[head]

example 3:
suppose u = [1, 2, 3]; v = [1, -1], then r = conv(u, v, 'valid') = [1, 1]
the procedures is as follows
1. flip v, so we have v = [-1, 1]
2. align u, v. however in each advancement of v, it must ensure u and v are totally overlapped. the precedures will automatically stop if partial overlap happen.

conv2: convolution of 2-D matrices U and V
suppose 
U = [1, 2, 3;
        4, 5, 6;
        7, 8, 9;];
V = [-1, 1;
        -1, 1];

R = conv2(U, V) = [ -1,  -1,  -1,  3;
                                -5,  -2,  -2,   9;
                               -11, -2,  -2, 15;
                                -7,  -1,  -1,   9;];

R = conv2(U, V, 'same') = [ -2,  -2,   9;
                                            -2,  -2, 15;
                                            -1,  -1,   9;];
 
R = conv2(U, V, 'valid') = [ -2,  -2;
                                           -2,  -2;];
the detailed procedures are identically the case of one dimension. note the way of reversing V

suppose V = [1, 2;
                      3, 4;];
then the reversing of V equals to
[4, 3;
 2, 1;]


filter2: filtering one two-dimensional vector using a FIR digital filter with coefficients stored in another two-dimensional vector

U = [1, 2, 3;
        4, 5, 6;
        7, 8, 9;];
V = [-1, 1;
        -1, 1];

R = filter2(V, U) = [ 2,  2,   -9;
                               -2,  2, -15;
                                1,  1,   -9;];
 
note: the filter2 function is essential resemble the conv2 function. and the filtering matrix needn't be reversed. 
 
however, pay attention to the order of input parameters. since the default behaviour is with shape parameter set as 'valid', so if U and V are of different sizes (usually they are of different sizes), different outputs will be generated if reverse the two.
 
gradient: find the gradient for 1-D vector or 2-D matrix
example 1: suppose a = [1, 2, 4, 1], g = gradient(a) = [1, 1.5, -0.5, -3];
the procedures are as follows:
suppose a[0] = 1, a[1] = 2, a[2] = 4, a[1] = 1;
so g[0] = (a[1] - a[0]) / 1 = 1; g[1] = (a[2] - a[0]) / 2 = 1.5;
g[2] = (a[3] - a[1]) / 2 = -0.5; g[3] = (a[3] - a[2]) / 1 = -3.
that's means for interior point, the central difference will be utilized. for edge point, the single-sided difference will be utilized.

example 2:
suppose a = [1, 2, 3;
                      4, 5, 6;
                      7, 8, 9;];
then [gx, gy] = gradient(a) leads to
gx = [1, 1, 1;
          1, 1, 1;
          1, 1, 1;];
gy = [3, 3, 3;
          3, 3, 3;
          3, 3, 3;];

pay attention to the ways of difference approaches applied to interior points and edge points along row and column directions respectively.

fspecial: create special 2-D filters, like gaussian, laplacian, etc.
example1: 
G = fspecial('gaussian', 5) = [
    0.0000    0.0000    0.0002    0.0000    0.0000
    0.0000    0.0113    0.0837    0.0113    0.0000
    0.0002    0.0837    0.6187    0.0837    0.0002
    0.0000    0.0113    0.0837    0.0113    0.0000
    0.0000    0.0000    0.0002    0.0000    0.0000
];
example2:
L = fspecial('laplacian') = [
    0.1667    0.6667    0.1667
    0.6667   -3.3333    0.6667
    0.1667    0.6667    0.1667
];
for other types, please refer to the help system of matlab
 
 
 
Reference:
UCF Computer Vision lectures
Signals and Systems for Bioengineers, John Semmlow
 
 
 
 
 

 








Sunday, September 20, 2015

Redirection of IO stream of C++

Sometimes, for programs written in C++, at the development stage, it's probably convenient to output the diagnosis information directly to console. And upon release, by the virtue of Macro, the diagnosis information can be eliminated. However for some daemon program, probably it still needs to output something if something wrong occurs. However simply turning on the Macro might not work since the TTY has been detached from standard output terminal or stand error terminal. The following are piece of work which can direct the IO stream from terminal to files. Hope it will be helpful under some occasion.


#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <string>
#include <sstream>
#include <iostream>

using namespace std;

class MyLog
{
public:
    MyLog(char* path)
    {
        m_fp = fopen(path, "w");
    }
    ~MyLog()
    {        
        if (m_fp)
        {
            fclose(m_fp);
        }
    }

    void operator() (const string& s)
    {
        if (m_fp)
        {
            fprintf(m_fp, "%s\n", s.c_str());
            fflush(m_fp);
        }
    }

    void operator() (const char* fmt, ...)
    {
        if (m_fp)
        {
            char buf[256];
            va_list ap;

            memset(buf, 0, sizeof(buf));

            va_start(ap, fmt);
            buf[vsnprintf(buf, sizeof(buf), fmt, ap)] = '\0';
            va_end(ap);

            fputs(buf, m_fp);
            fflush(m_fp);
        }
    }
    
private:
    FILE* m_fp;
};

MyLog myLog("/opt/test.txt");

class Redirect
{
public:
    Redirect(MyLog& log) : m_log(log), m_oldCoutStreamBuf(cout.rdbuf())
    {
        cout.rdbuf(m_strCout.rdbuf());
    }
    ~Redirect()
    {
        m_log(m_strCout.str());
        cout.rdbuf(m_oldCoutStreamBuf);    
    }

private:
    MyLog& m_log;
    streambuf* m_oldCoutStreamBuf;
    ostringstream m_strCout;
};

void func1()
{
    Redirect r(myLog);
    cout << "in func1!" << endl;
}

void func2()
{
    Redirect r(myLog);
    cout << "in func2!" << endl;
    func1();
}

int main(int argc, char* argv[])
{
    Redirect* p = new Redirect(myLog);
    cout << "in main!" << endl;
    func2();
    delete p;
    cout << "You should see this!" << endl;
    return 0;
}

  

Utilize Python to control GDB

There's a good article details how to control GDB via Python instead of GDB script, please refer: http://www.linuxjournal.com/article/11027.

However due to time evolves, the Python scripts seems to stop working. The following is a modified one to compatible with the example giving in the article.


#!/usr/local/bin/python
from collections import defaultdict
# A dictionary of mutex:owner
mutexOwners = {}
# A dictionary of blocking mutex:(threads..)
threadBlockers = defaultdict(list)
# Print the threads
#print "Process threads : "
#gdb.execute("info threads")
print "Analysing mutexes..."
# Step through processes running under gdb
def get_frame(hint):
    
    frame = gdb.selected_frame()

    if hint == None:
        return frame
    elif type(hint) == int:
        i = 0
        while i < hint and frame != None:
            pre_frame = frame.older()
            if frame == pre_frame and frame < hint - 1:
                frame = None
                break
            else:
                frame = new_frame
            i += 1
    elif type(hint) == str:
        while frame != None:
            pre_frame = frame.older()
            if frame.name() == hint or pre_frame == frame:
                break
            else:
                frame = pre_frame
    return frame

for process in gdb.inferiors():

    # Step through each thread in the process
    for thread in process.threads():

        # Examine the thread -- is it blocking on a mutex?
        thread.switch()
        frame = get_frame("pthread_mutex_lock")
        
        if frame != None:
            frame.select()

            # Make a note of which thread blocks on which mutex
            mutex = int(gdb.parse_and_eval("$r8"))
            threadBlockers[mutex].append(thread.ptid[0])

            # Make a note of which thread owns this mutex
            if not mutex in mutexOwners:
                v = gdb.parse_and_eval("*$r8").cast(gdb.lookup_type("pthread_mutex_t"))
                mutexOwners[mutex] = v['__data']['__owner']

    # Print the results of the analysis
    for mutex in mutexOwners:
        print "  Mutex 0x%x :" % mutex
        print "     -> held by thread : %d" % mutexOwners[mutex]
        s = ["%d" % t for t in threadBlockers[mutex]]
        print "     -> blocks threads : %s" % ' '.join(s)