Brief principle and VC Implementation of the musica (multi-scale image contrast enhancement) Algorithm

Musica patent documents: musica_patent-original document.pdf

Download the program source code: iplab_musica.rar

 

Algorithm principle:
The general methods of image enhancement are contrast stretching and histogram balancing. These two methods are suitable for the case where the gray level is too concentrated and a large number of gray levels are not fully utilized, the two methods are based on Histogram transformation and are independent of the pixel location information. If there is an original signal such as (a), we can see that the detailed signal (such as the high frequency, such as the fingerprint, clothes and wrinkles) is concentrated in a narrow gray range, and it is difficult for human eyes to distinguish. In addition, the intensity of detailed signal sets is within the entire gray-level space. Contrast enhancement and histogram balancing cannot be used to zoom in or out pixels in different regions to varying degrees. As a result, musica (multi-scale image contrast amplification) can be translated into multi-scale image contrast enhancement.
The basic principle is to extract the detailed information first, then enhance and enlarge the detailed signal, and then reconstruct it into the source image. The extraction of local details has always been a hot topic in research.
For the time being, no matter the theory, first from the perspective of practical use. (B) If smooth filtering is performed on the signal, a new blue signal can be obtained. It can be considered that the blue signal is similar to the low resolution of the original signal. It retains the overall features and discards the high-frequency detailed signal, for example, to analyze the overall features of an image, we need to use approximate signals to eliminate high-frequency (possibly noise) interference.
So (c), the original signal minus the approximate signal, the obtained is the high-frequency signal, that is, the detailed signal (of course, it may also be noise ).
Then, enhance (c). The simplest is to enlarge 2 times or perform log transformation to get (d ).
Then add (d) back to the approximate signal and reconstruct the original signal after enhancement, as shown in (e. We can see that the overall characteristics of the image have not changed, and the detailed signals of the two detail groups have been amplified, which is suitable for human eyes to recognize.

Figure 1

This is the most intuitive way to consider algorithm principles. It is just an approximate abstraction. The actual process requires more specific considerations. At the same time, this is my personal understanding, there will be some bias, and it remains to be improved.

Then, we have to give some specific principles and actual steps of the algorithm. The detailed description can be obtained from the musica patent document. Here is a personal understanding and description.
First, let's start with wavelet. intuitively, wavelet decomposition splits (w * h) images into four (W/2 * H/2) subgraphs, one of them is an approximate image, and the third is a wavelet coefficient matrix. Then, the approximate image is decomposed at the next level. For example, a is a low-resolution approximate image of the source image, D is a wavelet coefficient, and the physical meaning of D1, D2, and D3 can be viewed as a matrix of details in the horizontal, vertical, and diagonal directions. The reconstruction of wavelet is the inverse process of this decomposition process. We use four low-resolution subgraphs to restructure them into a high-resolution graph. It is best to look at the textbooks of wavelet to fully understand the characteristics of wavelet.

Figure 2

Intuitively looking at wavelet decomposition, the approximate image is a low-resolution approximation of the original image, retaining the overall features of the image, and the wavelet coefficient matrix is the difference between the low-resolution image and the original image, it is a high-frequency signal, that is, details. Through such decomposition, the details of each resolution can be obtained, so that features that are difficult to distinguish under a certain resolution can be easily divided and detected under another resolution. Human eyes can be automatically focused to select the best resolution to distinguish the current scene, and wavelet is also a lot of attention because it is similar to the human eye. Another important concept is wavelet decomposition. It can be understood as the wavelet basis used for image decomposition. Intuitively, the wavelet basis for different characteristics is similar to the same image, while the extracted detail matrix reflects the differences in image details. The selection of wavelet basis and the physical meanings of various wavelet features have not been well understood. One interpretation is that a certain wavelet base is representative of the wavelet matrix obtained from images of certain types of features. For example, XXX1 wavelet has a better effect on simple image combinations, xxx2 wavelet is good for noise-induced natural images and xxx3 is good for striped images... this is only an assumption discussed by one student and is to be verified.
However, this decomposition is not used in musica, but the Gaussian image pyramid is used in Rafael C. gonzarez. digital image processing. (Electronic Industry Press. 2004) I have introduced this book. Here I will refer to two of them to illustrate that the lower wave decomposition mentioned above serves as the foundation of mathematics. There are many teaching materials and it is easy to go deeper.
Gaussian image pyramid is also a kind of wavelet decomposition, which generates a low resolution approximation of W/2 * H/2, and a large residual array (that is, high-frequency details) such as the source image ). The sampling in the figure (2 H) refers to the use of the W * H image to obtain a subgraph of W/2 * H/2. The interpolation tool (2 rows) inserts a value between the row and column of the W/2 * H/2 image and expands it to the W * H image. The inserted value can be 0, or the mean of the left and right rows. Approximate filtering and insert filtering are the filters converted from wavelet basis. Gaussian filter is selected here. The filter parameters provide a set of empirical parameters in the patent documents, however, on a specific device, it is necessary to re-calculate it based on the noise of the actual imaging to achieve optimal results.

Figure 3

The overall algorithm diagram provided by the patent documentation is as follows:

Figure 4

2 is the input image, 30 is the decomposition process of the J-level image pyramid, and 31 is the residual matrix, and 31 is the J-level approximation of the image. 61 cache images (meaningful to specific devices) and 62 and perform the S-type function lookup transformation for 31. 34. In the reconstruction process, the transformed residual array 33 and J-level approximate 31 'are combined into the reinforced Image 4 and output.
Based on this process, you can start to consider algorithm program implementation.
Before programming on a PC, you must specify the following points:
The first is the effect of the algorithm. One is the choice of S-type functions in 62. Enhancement algorithms all have the side effects of noise enhancement. This problem is discussed in the patent documents, and the exponential transformation of Y = m * (x/m) P is selected, the value of P will affect the enhanced frequency band. Reasonable Selection of P and segmented selection of different P can suppress noise enhancement. The second is the 30 Decomposition Process and 34 reconstruction process. This is actually part of the wavelet decomposition. The series, decomposition steps, and selection of wavelet basis (here is the selection of filter parameters ), it will have an impact on the decomposition of the residual range array. Whether the residual range array has valid details is representative, and whether it can effectively filter noise is the key here. Some parameters can be set as variables in programming. The algorithm research and implementation can be separated by using the input in the preparation file. These two aspects are also discussed in the patent document for further reference.
It should be noted that it may not be appropriate to implement this algorithm on a PC. Generally, this algorithm is directly implemented on medical devices (embedded systems or some special firmware ), you can use the keyword "AGFA Musica" to find related information on Google. However, this is not taken into consideration here. We only use Musica as an interesting algorithm question and work together with some image processing knowledge to obtain the optimal solution. If you have any background, I hope you can introduce the relevant professional knowledge, which is more interesting.

Download the program source code: iplab_musica.rar

The overall process is relatively simple. The problem is some details in image processing, such as boundary processing during filtering and W/2 (H/2) When W or H is an odd number during image sampling) is the integer. There are also some problems with memory allocation during programming, and the allocation of temporary cache. However, if you do not consider the complexity, you can use the most common method. The above section only lists the approximate problems. Pay attention to them when encountering these problems.
Pseudo code of the entry function:
Void gdiplusimage: ipfuncmusica ()
{
O = original image;

For (decomposition level)
{Store the length and width of each level to the matrix s ;}

Imagegrey24to8 (o); // grayscale, only grayscale images are processed

X = musica_decomposition (O, S); // decomposition. X indicates that all images are arranged in the order of S, 30

X = musica_mapping (x, S); // transform the decomposed image, 32

O = musica_reconstruction (x, S); // refactor to obtain the enhanced image, 34

Imagegrey8to24 (o); // converts the image to 24, which is convenient for display and subsequent processing.
}

The musica_mapping (x, S) function is relatively simple. You only need to replace the gray-level look-up table for the residual array in S. The replaced table is constructed in advance based on the S-type function, and the process is similar to reverse color processing, refer to the previous article. You can consider separating the structure of the replace table and passing parameters to facilitate improvement.

Musica_decomposition decomposition and musica_reconstruction refactoring are relatively complex, and several improved decomposition and reconstruction processes are provided in the Patent Documents. Programming implementation is actually essentially the same. Therefore, the attached program implements the 4.a and 4.bin the patent document, which is the simplest one. Only the implementation of the 4.a decomposition is described, and the 4.c reconstruction is the inverse process. Refer to the source code.

Figure 5

One decomposition process:
1) Perform a first-level decomposition on X, and copy X to the buffer.
2) perform smooth filtering on the buffer.
3) The approximate graph X + 1 of the next level is obtained by buffer interlace sampling.
4) approximation graph X + 1 interlace interpolation 0 to buffer
5) perform smooth filtering by multiplying the filter parameter 4 (compensating for the row and column with 0)
6) X minus the buffer after smoothing to get the residual image stored in X

In this way, the next level of x + 1 is decomposed.

Figure 6

The program is stored in XP +. compile and run the program in net2003. After the program runs, open the image in the project directory and click shortcut 1 or the menu item "Image Processing-> Musica". Then, wait a moment and the processing speed is slow. Two sample images are shown in the project. The figure is the original image, and the figure is displayed as the result after the processing by a teacher on a professional device. Of course, the program running result is much worse than it, but the basic direction is correct. There are still many improvement steps and parameter merit-based in the patent documents, and interested friends can do it on their own. :)