Something More for Research

Explorer of Research #HEMBAD

Archive for the ‘Image / Video Filters’ Category

OpenCV 3.1 with CUDA , QT , Python Complete Installation on Windows in With Extra Modules

Posted by Hemprasad Y. Badgujar on May 13, 2016


OpenCV 3.1 with CUDA , QT , Python Complete Installation on Windows in With Extra Modules

The description here was tested on Windows 8.1 Pro. Nevertheless, it should also work on any other relatively modern version of Windows OS. If you encounter errors after following the steps described below, feel free to contact me.

Note :  To use the OpenCV library you have two options: Installation by Using the Pre-built Libraries or Installation by Making Your Own Libraries from the Source Files. While the first one is easier to complete, it only works if you are coding with the latest Microsoft Visual Studio IDE and doesn’t take advantage of the most advanced technologies we integrate into our library.

I am going to skip Installation by Using the Pre-built Libraries is it is easier to install even for New User. So Let’s Work on Installation by Making Your Own Libraries from the Source Files (If you are building your own libraries you can take the source files from OpenCV Git repository.) Building the OpenCV library from scratch requires a couple of tools installed beforehand

Prerequisites  Tools

Step By Step Prerequisites Setup

  • IDE : Microsoft Visual Studio. However, you can use any other IDE that has a valid CC++ compiler.

    Installing By Downloading from the Product Website Start installing Visual Studio by going to Visual Studio Downloads on the MSDN website and then choosing the edition you want to download.here we will going to use Visual Studio 2012 / ISO keys with  Visual Studio 2012 Update 4 /ISO and Step By Step Installing Visual Studio Professional 2012
  • Make Tool : Cmake is a cross-platform, open-source build system.

    Download and install the latest stable binary version: here we will going to use CMake 3 Choose the windows installer (cmake-x.y.z-win32.exe) and install it. Letting the cmake installer add itself to your path will make it easier but is not required.

  • Download OpenCV source Files by GIT/Sourceforge by : TortoiseGit or download source files from page on Sourceforge.

    The Open Source Computer Vision Library has >2500 algorithms, extensive documentation and sample code for real-time computer vision. It works on Windows, Linux, Mac OS X, Android and iOS.

  •  Python and Python libraries : Installation notes

    • It is recommended to uninstall any other Python distribution before installing Python(x,y)
    • You may update your Python(x,y) installation via individual package installers which are updated more frequently — see the plugins page
    • Please use the Issues page to request new features or report unknown bugs
    • Python(x,y) can be easily extended with other Python libraries because Python(x,y) is compatible with all Python modules installers: distutils installers (.exe), Python eggs (.egg), and all other NSIS (.exe) or MSI (.msi) setups which were built for Python 2.7 official distribution – see the plugins page for customizing options
    • Another Python(x,y) exclusive feature: all packages are optional (i.e. install only what you need)
    • Basemap users (data plotting on map projections): please see the AdditionalPlugins
  • Sphinx is a python documentation generator

    After installation, you better add the Python executable directories to the environment variable PATH in order to run Python and package commands such as sphinx-build easily from the Command Prompt.

    1. Right-click the “My Computer” icon and choose “Properties”

    2. Click the “Environment Variables” button under the “Advanced” tab

    3. If “Path” (or “PATH”) is already an entry in the “System variables” list, edit it. If it is not present, add a new variable called “PATH”.

      • Right-click the “My Computer” icon and choose “Properties”

      • Click the “Environment Variables” button under the “Advanced” tab

      • If “Path” (or “PATH”) is already an entry in the “System variables” list, edit it. If it is not present, add a new variable called “PATH”.

      • Add these paths, separating entries by ”;”:

        • C:\Python27 – this folder contains the main Python executable
        • C:\Python27\Scripts – this folder will contain executables added by Python packages installed with pip (see below)

        This is for Python 2.7. If you use another version of Python or installed to a non-default location, change the digits “27” accordingly.

      • Now run the Command Prompt. After command prompt window appear, type python and Enter. If the Python installation was successful, the installed Python version is printed, and you are greeted by the prompt>>> . TypeCtrl+Z and Enter to quit.

        Add these paths, separating entries by ”;”:

        • C:\Python27 – this folder contains the main Python executable
        • C:\Python27\Scripts – this folder will contain executables added by Python packages installed with easy_install (see below)

        This is for Python 2.7. If you use another version of Python or installed to a non-default location, change the digits “27” accordingly.

        After installation, you better add the Python executable directories to the environment variable PATH in order to run Python and package commands such as sphinx-build easily from the Command Prompt.

      • Install the pip command

      Python has a very useful pip command which can download and install 3rd-party libraries with a single command. This is provided by the Python Packaging Authority(PyPA): https://groups.google.com/forum/#!forum/pypa-dev

      To install pip, download https://bootstrap.pypa.io/get-pip.py and save it somewhere. After download, invoke the command prompt, go to the directory with get-pip.py and run this command:

      C:\> python get-pip.py
      

      Now pip command is installed. From there we can go to the Sphinx install.

      Note :pip has been contained in the Python official installation after version

        of Python-3.4.0 or Python-2.7.9.
  • Installing Sphinx with pip

    If you finished the installation of pip, type this line in the command prompt:

    C:\> pip install sphinx
    

    After installation, type sphinx-build -h on the command prompt. If everything worked fine, you will get a Sphinx version number and a list of options for this command.

    That it. Installation is over. Head to First Steps with Sphinx to make a Sphinx project.

    Now run the Command Prompt. After command prompt window appear, type python and Enter. If the Python installation was successful, the installed Python version is printed, and you are greeted by the prompt>>> . TypeCtrl+Z and Enter to quit.

  • Install the easy_install command

    Python has a very useful easy_install command which can download and install 3rd-party libraries with a single command. This is provided by the “setuptools” project: https://pypi.python.org/pypi/setuptools.

    To install setuptools, download https://bitbucket.org/pypa/setuptools/raw/bootstrap/ez_setup.py and save it somewhere. After download, invoke the command prompt, go to the directory with ez_setup.py and run this command:

    C:\> python ez_setup.py
    

    Now setuptools and its easy_install command is installed. From there we can go to the Sphinx install.

    Installing Sphinx with easy_install

    If you finished the installation of setuptools, type this line in the command prompt:
    C:\> easy_install sphinx
    

    After installation, type sphinx-build on the command prompt. If everything worked fine, you will get a Sphinx version number and a list of options for this command.

  • Numpy is a scientific computing package for Python. Required for the Python interface.

Try the (unofficial) binaries in this site: http://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy
You can get numpy 1.6.2 x64 with or without Intel MKL libs to Python 2.7

I suggest WinPython, a Python 2.7 distribution for Windows with both 32- and 64-bit versions.

It is also worth considering the Anaconda Python distribution. http://continuum.io/downloads

  • Numpy :Required for the Python interface abve Installation of python alosho included with Numpy and Scipy libraries

  • Intel © Threading Building Blocks (TBB) is used inside OpenCV for parallel code snippets. Download Here

    • Download TBB
      • Go to TBB download page to download the open source binary releases. I choose Commercial Aligned Release, because this has the most stable releases. I downloaded tbb43_20141204os, TBB 4.3 Update 3, specifically tbb43_20141204os for Windows. The release has the header files as well as the import library and DLL files prebuilt for Microsoft Visual C++ 8.0 and 9.0 on both x86(IA32) and x64(intel64). If you are aggressive and need the source code of TBB, you can try stable releases or development releases.
    • Install TBB
      • Extract the files in the zip file to a local directory, for example, C:\TBB. You should find tbb22_013oss under it. This is the installation directory, and doc, example, include etc should be directly under the installation folder.
      • Set a Windows environment variable TBB22_INSTALL_DIR to the above directory, e.g., C:\TBB\tbb22_013oss.
    • Develop with TBB
      • Add $(TBB22_INSTALL_DIR)\include to your C++ project’s additional include directories.
      • Add $(TBB22_INSTALL_DIR)\<arch>\<compiler>\lib (e.g., $(TBB22_INSTALL_DIR)\ia32\vc9\lib) to your project’s additional library directories.
      • Add to your project’s additional dependencies tbb.lib (Release) or tbb_debug.lib (Debug).
      • Write your C++ code to use TBB. See code below as an example.
    • Deploy with TBB
      • The TBB runtime is in TBB DLLs (tbb.dll/tbbmalloc.dll/tbbmalloc_proxy.dll for Release, tbb_debug.dll/tbbmalloc_debug.dll/tbbmalloc_proxy_debug.dll for Debug). They can be found in $(TBB22_INSTALL_DIR)\\\bin.
      • Your executable should have these DLLs in the same folder for execution.

    Intel © Integrated Performance Primitives (IPP) may be used to improve the performance of color conversion.(Paid)

    Intel Parallel Studio XE 2015 – Cluster Edition includes everything in the Professional edition (compilers, performance libraries, parallel models, performance profiler, threading design/prototyping, and memory & thread debugger). It adds a MPI cluster communications library, along with MPI error checking and tuning to design, build, debug and tune fast parallel code that includes MPI.

  • Eigen is a C++ template library for linear algebra.

     “install” Eigen?

    In order to use Eigen, you just need to download and extract Eigen‘s source code (see the wiki for download instructions). In fact, the header files in the Eigen subdirectory are the only files required to compile programs using Eigen. The header files are the same for all platforms. It is not necessary to use CMake or install anything.

     simple first program

    Here is a rather simple program to get you started.

    #include <iostream>
    #include <Eigen/Dense>
     int main()
    {
    MatrixXd m(2,2);
    m(0,0) = 3;
    m(1,0) = 2.5;
    m(0,1) = -1;
    m(1,1) = m(1,0) + m(0,1);
    std::cout << m << std::endl;
    }
  • Installing CUDA Development Tools

    The setup of CUDA development tools on a system running the appropriate version of Windows consists of a few simple steps:

    • Verify the system has a CUDA-capable GPU.
    • Download the NVIDIA CUDA Toolkit.
    • Install the NVIDIA CUDA Toolkit.
    • Test that the installed software runs correctly and communicates with the hardware.
    • CUDA Toolkit will allow you to use the power lying inside your GPU. we will going to use CUDA 7.5  Toolkit

      To verify that your GPU is CUDA-capable, open the Control Panel (StartControl > Panel) and double click on System. In the System Propertieswindow that opens, click the Hardware tab, then Device Manager. Expand the Display adapters entry. There you will find the vendor name and model of your graphics card. If it is an NVIDIA card that is listed inhttp://developer.nvidia.com/cuda-gpus, your GPU is CUDA-capable.

      The Release Notes for the CUDA Toolkit also contain a list of supported products.

       Download the NVIDIA CUDA Toolkit

      The NVIDIA CUDA Toolkit is available at http://developer.nvidia.com/cuda-downloads.

      Choose the platform you are using and download the NVIDIA CUDA Toolkit

      The CUDA Toolkit contains the CUDA driver and tools needed to create, build and run a CUDA application as well as libraries, header files, CUDA samples source code, and other resources.

      Download Verification

      The download can be verified by comparing the MD5 checksum posted at http://developer.nvidia.com/cuda-downloads/checksums with that of the downloaded file. If either of the checksums differ, the downloaded file is corrupt and needs to be downloaded again.

      To calculate the MD5 checksum of the downloaded file, follow the instructions at http://support.microsoft.com/kb/889768.

      Install the CUDA Software

      Before installing the toolkit, you should read the Release Notes, as they provide details on installation and software functionality.

      Note: The driver and toolkit must be installed for CUDA to function. If you have not installed a stand-alone driver, install the driver from the NVIDIA CUDA Toolkit.

      Graphical Installation

      Install the CUDA Software by executing the CUDA installer and following the on-screen prompts.

      Silent Installation

      Alternatively, the installer can be executed in silent mode by executing the package with the -s flag. Additional flags can be passed which will install specific subpackages instead of all packages. Allowed subpackage names are: CUDAToolkit_6.5, CUDASamples_6.5, CUDAVisualStudioIntegration_6.5, and Display.Driver. For example, to install only the driver and the toolkit components:

      .exe -s CUDAToolkit_6.5 Display.Driver

      This will drastically improve performance for some algorithms (e.g the HOG descriptor). Getting more and more of our algorithms to work on the GPUs is a constant effort of the OpenCV .

  • JRE :Java run time environment

    Installing Ant The binary distribution of Ant consists of the following directory layout:

      ant
       +--- README, LICENSE, fetch.xml, other text files. //basic information
       +--- bin  // contains launcher scripts
       |
       +--- lib  // contains Ant jars plus necessary dependencies
       |
       +--- docs // contains documentation
       |      |
       |      +--- images  // various logos for html documentation
       |      |
       |      +--- manual  // Ant documentation (a must read ;-)
       |
       +--- etc // contains xsl goodies to:
                //   - create an enhanced report from xml output of various tasks.
                //   - migrate your build files and get rid of 'deprecated' warning
                //   - ... and more ;-)
    

    Only the bin and lib directories are required to run Ant. To install Ant, choose a directory and copy the distribution files there. This directory will be known as ANT_HOME.

Before you can run Ant there is some additional set up you will need to do unless you are installing the RPM version from jpackage.org:

  • Add the bin directory to your path.
  • Set the ANT_HOME environment variable to the directory where you installed Ant. On some operating systems, Ant’s startup scripts can guess ANT_HOME(Unix dialects and Windows NT/2000), but it is better to not rely on this behavior.
  • Optionally, set the JAVA_HOME environment variable (see the Advanced section below). This should be set to the directory where your JDK is installed.

Operating System-specific instructions for doing this from the command line are in the Windows, Linux/Unix (bash), and Linux/Unix (csh) sections. Note that using this method, the settings will only be valid for the command line session you run them in. Note: Do not install Ant’s ant.jar file into the lib/ext directory of the JDK/JRE. Ant is an application, whilst the extension directory is intended for JDK extensions. In particular there are security restrictions on the classes which may be loaded by an extension.

Windows Note:
The ant.bat script makes use of three environment variables – ANT_HOME, CLASSPATH and JAVA_HOME. Ensure that ANT_HOME and JAVA_HOME variables are set, and that they do not have quotes (either ‘ or “) and they do not end with \ or with /. CLASSPATH should be unset or empty.

Check Installation

You can check the basic installation with opening a new shell and typing ant. You should get a message like this

Buildfile: build.xml does not exist!
Build failed

So Ant works. This message is there because you need to write an individual buildfile for your project. With a ant -version you should get an output like

Apache Ant(TM) version 1.9.2 compiled on July 8 2013

If this does not work ensure your environment variables are set right. They must resolve to:

  • required: %ANT_HOME%\bin\ant.bat
  • optional: %JAVA_HOME%\bin\java.exe
  • required: %PATH%=…maybe-other-entries…;%ANT_HOME%\bin;…maybe-other-entries

ANT_HOME is used by the launcher script for finding the libraries. JAVA_HOME is used by the launcher for finding the JDK/JRE to use. (JDK is recommended as some tasks require the java tools.) If not set, the launcher tries to find one via the %PATH% environment variable. PATH is set for user convenience. With that set you can just start ant instead of always typingthe/complete/path/to/your/ant/installation/bin/ant.

Advertisements

Posted in GPU (CUDA), Image / Video Filters, Mixed, OpenCV, OpenCV | Leave a Comment »

Bilateral Filtering

Posted by Hemprasad Y. Badgujar on September 14, 2015


Popular Filters

When smoothing or blurring images (the most popular goal of smoothing is to reduce noise), we can use diverse linear filters, because linear filters are easy to achieve, and are kind of fast, the most used ones are Homogeneous filter, Gaussian filter, Median filter, et al.

When performing a linear filter, we do nothing but output pixel’s value g(i,j)  which is determined as a weighted sum of input pixel values f(i+k, j+l):

g(i, j)=SUM[f(i+k, j+l) h(k, l)];

in which, h(k, l)) is called the kernel, which is nothing more than the coefficients of the filter.

Homogeneous filter is the most simple filter, each output pixel is the mean of its kernel neighbors ( all of them contribute with equal weights), and its kernel K looks like:

1

 Gaussian filter is nothing but using different-weight-kernel, in both x and y direction, pixels located in the middle would have bigger weight, and the weights decrease with distance from the neighborhood center, so pixels located on sides have smaller weight, its kernel K is something like (when kernel is 5*5):

gkernel

Median filter is something that replace each pixel’s value with the median of its neighboring pixels. This method is great when dealing with “salt and pepper noise“.

Bilateral Filter

By using all the three above filters to smooth image, we not only dissolve noise, but also smooth edges, which make edges less sharper, even disappear. To solve this problem, we can use a filter called bilateral filter, which is an advanced version of Gaussian filter, it introduces another weight that represents how two pixels can be close (or similar) to one another in value, and by considering both weights in image,  Bilateral filter can keep edges sharp while blurring image.

Let me show you the process by using this image which have sharp edge.

21

 

Say we are smoothing this image (we can see noise in the image), and now we are dealing with the pixel at middle of the blue rect.

22   23

Left-above picture is a Gaussian kernel, and right-above picture is Bilateral filter kernel, which considered both weight.

We can also see the difference between Gaussian filter and Bilateral filter by these pictures:

Say we have an original image with noise like this

32

 

By using Gaussian filter, the image is smoother than before, but we can see the edge is no longer sharp, a slope appeared between white and black pixels.

33

 

However, by using Bilateral filter, the image is smoother, the edge is sharp, as well.

31

OpenCV code

It is super easy to make these kind of filters in OpenCV:

1 //Homogeneous blur:
2 blur(image, dstHomo, Size(kernel_length, kernel_length), Point(-1,-1));
3 //Gaussian blur:
4 GaussianBlur(image, dstGaus, Size(kernel_length, kernel_length), 0, 0);
5 //Median blur:
6 medianBlur(image, dstMed, kernel_length);
7 //Bilateral blur:
8 bilateralFilter(image, dstBila, kernel_length, kernel_length*2, kernel_length/2);

and for each function, you can find more details in OpenCV Documentation

Test Images

Glad to use my favorite Van Gogh image :

vangogh

 

From left to right: Homogeneous blur, Gaussian blur, Median blur, Bilateral blur.

(click iamge to view full size version :p )

kernel length = 3:

homo3 Gaussian3 Median3 Bilateral3

kernel length = 9:

homo9 Gaussian9 Median9 Bilateral9
kernel length = 15:

homo15 Gaussian15 Median15 Bilateral15

kernel length = 23:

homo23 Gaussian23 Median23 Bilateral23
kernel length = 31:

homo31 Gaussian31 Median31 Bilateral31
kernel length = 49:

homo49 Gaussian49 Median49 Bilateral49
kernel length = 99:

homo99 Gaussian99 Median99 Bilateral99

Trackback URL.

Posted in C, Image / Video Filters, Image Processing, OpenCV, OpenCV, OpenCV Tutorial | Leave a Comment »

Display more than one image in a single window

Posted by Hemprasad Y. Badgujar on February 6, 2015


There is no inbuilt support to display more than one image in OpenCV. Here is a function illustrating how to display more than one image in a single window using Intel OpenCV. The method used is to set the ROIs of a Single Big image and then resizing and copying the input images on to the Single Big Image.

  1#include <cv.h>
  2#include <highgui.h>
  3
  4#include <stdio.h>
  5#include <stdarg.h>
  6
  7/*Function///////////////////////////////////////////////////////////////
  8
  9Name:       cvShowManyImages
 10
 11Purpose:    This is a function illustrating how to display more than one 
 12               image in a single window using Intel OpenCV
 13
 14Parameters: char *title: Title of the window to be displayed
 15            int nArgs:   Number of images to be displayed
 16            ...:         IplImage*, which contains the images
 17
 18Language:   C++
 19
 20The method used is to set the ROIs of a Single Big image and then resizing 
 21and copying the input images on to the Single Big Image.
 22
 23This function does not stretch the image... 
 24It resizes the image without modifying the width/height ratio..
 25
 26This function can be called like this:
 27
 28cvShowManyImages("Images", 2, img1, img2);
 29or
 30cvShowManyImages("Images", 5, img2, img2, img3, img4, img5);
 31
 32This function can display upto 12 images in a single window.
 33It does not check whether the arguments are of type IplImage* or not.
 34The maximum window size is 700 by 660 pixels.
 35Does not display anything if the number of arguments is less than
 36    one or greater than 12.
 37
 38If you pass a pointer that is not IplImage*, Error will occur.
 39Take care of the number of arguments you pass, and the type of arguments, 
 40which should be of type IplImage* ONLY.
 41
 42Idea was from [[BettySanchi]] of OpenCV Yahoo! Groups.
 43
 44If you have trouble compiling and/or executing
 45this code, I would like to hear about it.
 46
 47You could try posting on the OpenCV Yahoo! Groups
 48[url]http://groups.yahoo.com/group/OpenCV/messages/ [/url]
 49
 50Parameswaran, 
 51Chennai, India.
 52
 53cegparamesh[at]gmail[dot]com            
 54
 55...
 56///////////////////////////////////////////////////////////////////////*/
 57
 58void cvShowManyImages(char* title, int nArgs, ...) {
 59
 60    // img - Used for getting the arguments 
 61    IplImage *img;
 62
 63    // [[DispImage]] - the image in which input images are to be copied
 64    IplImage *DispImage;
 65
 66    int size;
 67    int i;
 68    int m, n;
 69    int x, y;
 70
 71    // w - Maximum number of images in a row 
 72    // h - Maximum number of images in a column 
 73    int w, h;
 74
 75    // scale - How much we have to resize the image
 76    float scale;
 77    int max;
 78
 79    // If the number of arguments is lesser than 0 or greater than 12
 80    // return without displaying 
 81    if(nArgs <= 0) {
 82        printf("Number of arguments too small....\n");
 83        return;
 84    }
 85    else if(nArgs > 12) {
 86        printf("Number of arguments too large....\n");
 87        return;
 88    }
 89    // Determine the size of the image, 
 90    // and the number of rows/cols 
 91    // from number of arguments 
 92    else if (nArgs == 1) {
 93        w = h = 1;
 94        size = 300;
 95    }
 96    else if (nArgs == 2) {
 97        w = 2; h = 1;
 98        size = 300;
 99    }
100    else if (nArgs == 3 || nArgs == 4) {
101        w = 2; h = 2;
102        size = 300;
103    }
104    else if (nArgs == 5 || nArgs == 6) {
105        w = 3; h = 2;
106        size = 200;
107    }
108    else if (nArgs == 7 || nArgs == 8) {
109        w = 4; h = 2;
110        size = 200;
111    }
112    else {
113        w = 4; h = 3;
114        size = 150;
115    }
116
117    // Create a new 3 channel image
118    [[DispImage]] = cvCreateImage( cvSize(100 + size*w, 60 + size*h), 8, 3 );
119
120    // Used to get the arguments passed
121    va_list args;
122    va_start(args, nArgs);
123
124    // Loop for nArgs number of arguments
125    for (i = 0, m = 20, n = 20; i < nArgs; i++, m += (20 + size)) {
126
127        // Get the Pointer to the IplImage
128        img = va_arg(args, IplImage*);
129
130        // Check whether it is NULL or not
131        // If it is NULL, release the image, and return
132        if(img == 0) {
133            printf("Invalid arguments");
134            cvReleaseImage(&DispImage);
135            return;
136        }
137
138        // Find the width and height of the image
139        x = img->width;
140        y = img->height;
141
142        // Find whether height or width is greater in order to resize the image
143        max = (x > y)? x: y;
144
145        // Find the scaling factor to resize the image
146        scale = (float) ( (float) max / size );
147
148        // Used to Align the images
149        if( i % w == 0 && m!= 20) {
150            m = 20;
151            n+= 20 + size;
152        }
153
154        // Set the image ROI to display the current image
155        cvSetImageROI(DispImage, cvRect(m, n, (int)( x/scale ), (int)( y/scale )));
156
157        // Resize the input image and copy the it to the Single Big Image
158        cvResize(img, DispImage);
159
160        // Reset the ROI in order to display the next image
161        cvResetImageROI(DispImage);
162    }
163
164    // Create a new window, and show the Single Big Image
165    cvNamedWindow( title, 1 );
166    cvShowImage( title, DispImage);
167
168    cvWaitKey();
169    cvDestroyWindow(title);
170
171    // End the number of arguments
172    va_end(args);
173
174    // Release the Image Memory
175    cvReleaseImage(&DispImage);
176}
177

You can use this function as in this sample program:

 1int main() {
 2
 3    IplImage *img1 = cvLoadImage("Image01.jpg");
 4
 5    IplImage *img2 = cvLoadImage("Image02.jpg");
 6
 7    IplImage *img3 = cvLoadImage("Image03.jpg");
 8
 9    IplImage *img4 = cvLoadImage("Image04.jpg");
10
11    IplImage *img5 = cvLoadImage("Image05.jpg");
12
13    IplImage *img6 = cvLoadImage("Image06.jpg");
14
15    cvShowManyImages("Image", 6, img1, img2, img3, img4, img5, img6);
16
17    return 0;
18
19}
20

The method used is to set the ROIs of a Single Big image and then resizing and copying the input images on to the Single Big Image.

This function does not stretch the image… It resizes the image without modifying the width/height ratio..

This function can be called like this:

1   cvShowManyImages("Image", 2, img1, img2);

or

1   cvShowManyImages("Image", 6, img1, img2, img3, img4, img5, img6);

upto 12 images.

This function can display upto 12 images in a single window. It does not check whether the arguments are of type IplImage* or not. The maximum window size is 700 by 660 pixels. Does not display anything if the number of arguments is less than one or greater than 12.

If you pass a pointer that is not IplImage*, Error will occur. Take care of the number of arguments you pass, and the type of arguments, which should be of type IplImage* ONLY.

Idea was from BettySanchi of OpenCV Yahoo! Groups.

If you have trouble compiling and/or executing this code, I would like to hear about it.

Here is a sample ScreenShot:

Posted in Computer Vision, GPU (CUDA), Image / Video Filters, OpenCV, PARALLEL | Tagged: , | Leave a Comment »

How to run CUDA 6.5 in Emulation Mode

Posted by Hemprasad Y. Badgujar on December 20, 2014


How to run CUDA in Emulation Mode

Some beginners feel a little bit dejected when they find that their systems donotcontainGPUs to learn andworkwithCUDA. In this blog post, I shall include the step by step process of installingandexecutingCUDA programs in emulation mode on a system with no GPU installed in it.It is mentioned here thatyouwill not be able to gain any performance advantage expected out of a GPU (obviously). Instead, the performance will be worse than a CPU implementation. However, emulation mode provides an excellent tool to compile and debugyourCUDA codes for more advanced purposes.Please note that I performed the following steps for a Dell Xeon with Windows 7 (32-bit) system.1. Acquire and install Microsoft Visual Studio 2008 on your system.

2. Access the CUDA Toolkit Archives  page and select CUDA Toolkit 6.0 / 6.5 version. (It is the last version that came with emulation mode. Emulation mode was discontinued in later versions.)

3. Download and install the following on your machine:-

  • Developer Drivers for Win8/win7 X64  – (Select the one as required for your machine.)
  • CUDA Toolkit
  • CUDA SDK Code Samples
  • CUBLAS and CUFFT (If required)

4. The next step is to check whether the sample codes run properly on the system or not. This will ensure that there is nothing missing from the required installations. Browse the nVIDIA GPU Computing SDK using the windows start bar or by using the following path in your My Computer address bar:-
As per your working Platform
“C:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK\C\bin\win32\Release”
“C:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK\C\bin\win64\Release”

(Also note that the ProgramData folder is by default set to “Hidden” attribute. It will be good if you unhide theis folder as it will be frequently utilized later on as you progress with your CUDA learning spells.)

5. Run the “deviceQuery” program and it should output something similar as shown in Fig. 1. Upon visual inspection of the output data, it can be seen that “there is no GPU device found” however the test has PASSED. This means that all the required installations for CUDA in emulation mode has been completed and now we can proceed with writing, compiling and executing CUDA programs in emulation mode.

Figure 1. Successful Rxecution of deviceQuery.exe ** Demo Example only

6. Open Visual Studio and create a new Win32 console project. Let’s name it “HelloCUDAEmuWorld”. Remember to select the “EMPTY PROJECT” option in Application Settings. Now Right Click on “Source Files” in the project tree and add new C++ code item. Remember to include the extension “.cu” instead of “.cpp”. Let’s name this item as “HelloCUDAEmuWorld.cu”. (If you forget the file extension, it can always be renamed via the project tree on the left).

7. Include the CUDA include, lib and  bin paths to MS Visual Studio. They were located at “C:\CUDA” in my system.

The next steps need to be performed for every new CUDA project when created.

8. Right Click on the project and select Custom Build Rules. Check the Custom Build Rules v6.0.0 option if available. Otherwise, click on Find Existing… and navigate to “C:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK\C\common” and select Cuda.rules. This will add the build rules for CUDA v6.0to VS 2012.

9. Right click on the project and select Properties. Navigate to Configuration Properties –> Linker –> Input. Type in cudart.lib in the Additional Dependencies text bar and click Okay. Now we are ready to compile and run our first ever CUDA program in emulation mode. But first we need to activate the emulation  mode for .cu files.

10. Once again  Right click on the project and select Properties. Navigate to Configuration Properties –> CUDA Build Rule v6.0.0 –> General. Set Emulation Mode from No to Yes in the right hand column of the opened window. Click Okay.

11. Type in the following in the code editor and build and compile the project. And there it is. Your first ever CUDA program, in Emulation Mode. Something to brag about among friends.

int main(void)
{
return 0;
}

I hope this effort would not go in vain and offer some help to anyone whois tied upregarding this issue. Do contact if there is any queryregarding the above procedure.Source (http://hpcden.blogspot.in)

Posted in Computer Vision, Computing Technology, CUDA, GPU (CUDA), GPU Accelareted, Image / Video Filters, My Research Related, OpenCV, PARALLEL, Project Related | Leave a Comment »

Posted by Hemprasad Y. Badgujar on October 12, 2014


Mathematics Function handbook

Abstract. Function handbook contains all elementary, trigonometric and hyperbolic functions — their definitions, graphs, properties and identities.

 

Function handbook

1. Introduction

When working at scientific calculator Li-L project we wanted to give some gift to its the future users. This scientific calculator has an embedded Help available just in two clicks for every part of the application and we had an idea to create function handbook and include it as a part of the Help, so that click in the Li-L calculator keypad took the user to the handbook page with function description.

The most challenging part of the idea was graphics. We wanted every function to have its graph and we set high standards for graphics: it was to be properly scaled, numbered and antialiased. It was interesting and difficult job: every figure was created in several stages with final stage of editing on pixel level. When the task was completed, maybe first time in my life I saw all the functions in the same scale, properly curved and cut, so that all the features were clearly visible.

The result got was so great we decided to put this beauty online. And now you can access what Li-L calculator users have on their desktops. Below you will find functions put into several lists — every list holds all functions, but ordering is different, so that you can look up the function by its notation, name or article head.

Li-L calculator handbook contains all elementary, trigonometric and hyperbolic functions. Scientific calculator Li-X got handbook extended with gamma and Bessel functions and scientific calculator Li-Xc got handbook extended with modified Bessel functions.

2. Content

  1. or sqrt — square root
  2. ^ — power
  3. Γ — gamma function
  4. abs — absolute value
  5. arccos — trigonometric arc cosine
  6. arccot or arcctg — trigonometric arc cotangent
  7. arccsc or arccosec — trigonometric arc cosecant
  8. arcosh or arch — arc-hyperbolic cosine
  9. arcoth or arcth — arc-hyperbolic cotangent
  10. arcsch — arc-hyperbolic cosecant
  11. arcsec — trigonometric arc secant
  12. arcsin — trigonometric arc sine
  13. arctan or arctg — trigonometric arc tangent
  14. arsech or arsch — arc-hyperbolic secant
  15. arsinh or arsh — arc-hyperbolic sine
  16. artanh or arth — arc-hyperbolic tangent
  17. ceil — ceiling
  18. cos — trigonometric cosine
  19. cosh or ch — hyperbolic cosine
  20. cot or ctg — trigonometric cotangent
  21. coth or cth — hyperbolic cotangent
  22. csc or cosec — trigonometric cosecant
  23. csch — hyperbolic cosecant
  24. exp — exponent
  25. floor — floor
  26. Iν — modified Bessel function of the first kind
  27. Jν — Bessel function of the first kind
  28. Kν — modified Bessel function of the second kind
  29. ln — natural logarithm
  30. log or lg — decimal logarithm
  31. n! — factorial
  32. sec — trigonometric secant
  33. sech — hyperbolic secant
  34. sign — signum
  35. sin — trigonometric sine
  36. sinh or sh — hyperbolic sine
  37. tan or tg — trigonometric tangent
  38. tanh or th — hyperbolic tangent
  39. Yν — Bessel function of the second kind

3. Function index by name

  1. Absolute value abs
  2. Arc-hyperbolic cosecant arcsch
  3. Arc-hyperbolic cosine arcosh or arch
  4. Arc-hyperbolic cotangent arcoth or arcth
  5. Arc-hyperbolic secant arsech or arsch
  6. Arc-hyperbolic sine arsinh or arsh
  7. Arc-hyperbolic tangent artanh or arth
  8. Bessel function of the first kind Jν
  9. Bessel function of the second kind Yν
  10. Ceiling function ceil
  11. Decimal logarithm log or lg
  12. Exponent exp
  13. Factorial n!
  14. Floor function floor
  15. Gamma function Γ
  16. Hyperbolic cosecant csch
  17. Hyperbolic cosine cosh or ch
  18. Hyperbolic cotangent coth or cth
  19. Hyperbolic secant sech
  20. Hyperbolic sine sinh or sh
  21. Hyperbolic tangent tanh or th
  22. Modified Bessel function of the first kind Iν
  23. Modified Bessel function of the second kind Kν
  24. Natural logarithm ln
  25. Power ^
  26. Signum sign
  27. Square root or sqrt
  28. Trigonometric arc cosecant arccsc or arccosec
  29. Trigonometric arc cosine arccos
  30. Trigonometric arc cotangent arccot or arcctg
  31. Trigonometric arc secant arcsec
  32. Trigonometric arc sine arcsin
  33. Trigonometric arc tangent arctan or arctg
  34. Trigonometric cosecant csc or cosec
  35. Trigonometric cosine cos
  36. Trigonometric cotangent cot or ctg
  37. Trigonometric secant sec
  38. Trigonometric sine sin
  39. Trigonometric tangent tan or tg

4. Function index by notation

  1. — square root
  2. ^ — power
  3. Γ — gamma function
  4. abs — absolute value
  5. arccos — trigonometric arc cosine
  6. arccosec — trigonometric arc cosecant
  7. arccot — trigonometric arc cotangent
  8. arccsc — trigonometric arc cosecant
  9. arcctg — trigonometric arc cotangent
  10. arch — arc-hyperbolic cosine
  11. arcosh — arc-hyperbolic cosine
  12. arcoth — arc-hyperbolic cotangent
  13. arcsch — arc-hyperbolic cosecant
  14. arcsec — trigonometric arc secant
  15. arcsin — trigonometric arc sine
  16. arctan — trigonometric arc tangent
  17. arctg — trigonometric arc tangent
  18. arcth — arc-hyperbolic cotangent
  19. arsch — arc-hyperbolic secant
  20. arsech — arc-hyperbolic secant
  21. arsh — arc-hyperbolic sine
  22. arsinh — arc-hyperbolic sine
  23. artanh — arc-hyperbolic tangent
  24. arth — arc-hyperbolic tangent
  25. ceil — ceiling
  26. ch — hyperbolic cosine
  27. cos — trigonometric cosine
  28. cosec — trigonometric cosecant
  29. cosh — hyperbolic cosine
  30. cot — trigonometric cotangent
  31. coth — hyperbolic cotangent
  32. csc — trigonometric cosecant
  33. csch — hyperbolic cosecant
  34. ctg — trigonometric cotangent
  35. exp — exponent
  36. floor — floor
  37. Iν — modified Bessel function of the first kind
  38. Jν — Bessel function of the first kind
  39. Kν — modified Bessel function of the second kind
  40. lg — decimal logarithm
  41. ln — natural logarithm
  42. log — decimal logarithm
  43. n! — factorial
  44. sec — trigonometric secant
  45. sech — hyperbolic secant
  46. sh — hyperbolic sine
  47. sign — signum
  48. sin — trigonometric sine
  49. sinh — hyperbolic sine
  50. sqrt — square root
  51. tan — trigonometric tangent
  52. tanh — hyperbolic tangent
  53. tg — trigonometric tangent
  54. th — hyperbolic tangent
  55. Yν — Bessel function of the second kind

References. The handbook is a part of scientific calculator Li-L and scientific calculator Li-X products.

Posted in Apps Development, Computer Languages, Computer Research, Free Tools, Image / Video Filters, My Research Related, PARALLEL | Leave a Comment »

Fast Fourier transform — FFT

Posted by Hemprasad Y. Badgujar on October 12, 2014


Fast Fourier transform — FFT

Category. Digital signal processing (DSP) software development.

Abstract. The article is a practical tutorial for fast Fourier transform — FFT — understanding and implementation. Article contains theory, C++ source code and programming instructions. Popular Cooley-Tukey technique is considered.

References. Case study: ECG processing — R-peaks detection and 4D ultrasound volume image processing.

Download fast Fourier transform — FFT — C++ source code (zip, 4 Kb)

Fourier transform

1. Introduction to fast Fourier transform

Fast Fourier transformFFT — is speed-up technique for calculating discrete Fourier transformDFT, which in turn is discrete version of continuous Fourier transform, which indeed is origin for all its versions. So, historically continuous form of the transform was discovered, then discrete form was created for sampled signals and then algorithm for fast calculation of discrete version was invented.

2. Understanding FFT

First of all let us have a look at what Fourier transform is. Fourier transform is an integral of the form:

Fourier transform (1)

The transform operates in complex domain. Recall, that imaginary exponent could be written as:

Imaginary exponent (2)

For sampled function continuous transform (1) turns into discrete one:

Discrete Fourier transform (DFT) (3)

Expression (3) is discrete Fourier transformDFT. Here {f0, f1, … , fN-1} is input discrete function and {F0, F1, … , FN-1} is result of Fourier transform.

It is easily could be seen that to program DFT it is enough to write double loop and just calculate sums of products of input samples and imaginary exponents. The number of operations required is obviously of O(N2) order. But due to transform properties it is possible to reduce the number of operations to the order of O(Nlog2N). Now, let us explore tricks we can use to speed-up calculations.

Let us put N=8 and write down our DFT:

Discrete Fourier transform for N=8 (4)

Easily could be seen we can split the sum into two similar sums separating odd and even terms and factoring out the latter sum:

Transform separation (5)

Now we can split the sums in brackets again:

Transform separation (6)

Thus, we have 3 — log28 — levels of summation. The deepest one in parenthesis, the middle one in brackets and the outer one. For every level exponent multiplier for the second term is the same.

Transform separation (7)

And now the most important observation one can make to get speed-up: periodicity of the exponent multipliers.

Imaginary exponent periodicity (8)

For the exponent multiplier in parenthesis period is n=2, which means sums in parenthesis are exactly the same for n=0, 2, 4, 6 and for n=1, 3, 5, 7. It means on deepest level in parenthesis we need 4×2=8 — number of sums times period — sums in total. And note, since n=1, 3, 5, 7 corresponds to the half of the period π, exponent multiplier is the same as for n=0, 2, 4, 6 but with the opposite sign.

Imaginary exponent periodicity (9)

Indeed they are 1 for n=0, 2, 4, 6 and -1 for n=1, 3, 5, 7.

For the exponent multiplier in brackets period is n=4, which means we have the same sums for pairs n=0, 4; n=1, 5; n=2, 6 and n=3, 7. It means on the middle level in brackets we have 2×4=8 sums and the second half of them could be received again by changing sign in first half of them — due to the fact distance between n and n+2 is π. Thus, for n=0, 4 the factor is 1 and for n=2, 6 it is -1; for n=1, 5 it equals –i and forn=3, 7 it is i.

On the outer level we have just one sum for every transform component, and the period of the exponent multiplier is 8. Which gives us 1×8=8 sums and second half of them could be received by changing sign in the first half.

So, on every calculation level we have 8 sums. In terms of N it means we have log2N levels and N sums on every level, which gives usO(Nlog2N) order of number of operations. On the other hand having the constant number of sums on every level means we can process data in-place.

In summary, we have got fast Fourier transformFFT. Now it is time to develop step-by-step instruction list to be carved in code.

3. FFT algorithm

Having understanding of what features of DFT we are going to exploit to speed-up calculation we can write down the following algorithm:

  1. Prepare input data for summation — put them into convenient order;
  2. For every summation level:
    1. For every exponent factor of the half-period:
      1. Calculate factor;
      2. For every sum of this factor:
        1. Calculate product of the factor and the second term of the sum;
        2. Calculate sum;
        3. Calculate difference.

First step of reordering is putting input data into their natural order in (7): {f0, f1, f2, f3, f4, f5, f6, f7}→{f0, f4, f2, f6, f1, f5, f3, f7}. Since this new order was received as result of successive splitting terms into even and odd ones, in binary domain it looks like ordering based on bit greatness starting from the least significant bit — see table 1.

0 000 000 000 0
1 001 010 100 4
2 010 100 010 2
3 011 110 110 6
4 100 001 001 1
5 101 011 101 5
6 110 101 011 3
7 111 111 111 7

Table. 1. Reordering in binary domain.

This leads to “mirrored arithmetics” — result binary column in the mirror looks like {0, 1, 2, 3, 4, 5, 6, 7} — our initial order indeed. Thus, to reorder input elements it is enough just to count in mirrored manner. Since double mirroring gives again the same number, reordering reduces to swaps.

Summation levels include our parenthesis, brackets and outer level. In general case this leads to iterations on pairs, quadruples, octads and so on.

Further, we iterate on components of half-period, second half-period we are getting as result of taking differences instead of sums for the first half. Thus, for the deepest level of parenthesis period is 2 and half-period is 1, which means this cycle will be executed only once. For the second level period is 4 and half-period is 2 and cycle will be executed 2 times. In general case we have succession 1, 2, 4, 8, … for this cycle.

Factor calculation is calculation of our imaginary exponent. To restrict the number of trigonometric function calls (2) we use the recurrence:

Trigonometric recurrence (10)

Which is indeed expression

Trigonometric recurrence (11)

written in a tricky way not to lose significance for small δ — for this case cosδ≈1 and sinδδ, which tells us that for cosδ memory will be just packed with .999999(9) but for sinδ there will be much more useful information. Thus, (10) is just the way to eliminate cosδ from calculations. If you look back at (7) you will find, that for N=8 δ=π/2, π/4 — not a small numbers. But for transforms of much bigger Nδ=π/2, π/4, π/8, … up to 2π/N, for sure, could be very small.

The innermost loop looks for sums, where calculated imaginary exponent present, calculates product and takes sum and difference, which is the sum for the second half-period at π distance, where our exponent changes its sign but not the absolute value according to (9). To perform in-place processing we utilize the following scheme:

Fig. 1. Butterfly.Fig. 1. Butterfly.

This operation is elementary brick of in-place FFT calculation and usually is called butterfly. The bottom term is multiplied by imaginary exponent and then sum of the terms is stored in place of the upper term and difference is stored in place of the bottom term. General butterfly picture is depicted below — fig. 2.

Fig. 2. All butterflies for N=8.Fig. 2. All butterflies for N=8.

Elegant scheme. It is time to engrave this beauty in code. But before we delve into programming let us make a small digression: it is known thing that Fourier transform is a transfer to frequency domain, so, let us see how to be back from the realm of frequencies.

4. Inverse Fourier transform

Expression for inverse Fourier transform is

Inverse Fourier transform (12)

and its discrete counterpart is

Inverse discrete Fourier transform (13)

Thus, the difference between forward (3) and inverse (13) transforms is just a sign and not necessary scale factor — one does not need it if interested in ratio between components but not in absolute values. It means that the routine for forward transform with slight modification can perform inverse one as well.

5. FFT programming

Here is our FFT class declaration.

class CFFT
{
public:
   //   FORWARD FOURIER TRANSFORM
   //     Input  - input data
   //     Output - transform result
   //     N      - length of both input data and result
   static bool Forward(const complex *const Input, complex *const Output,
      const unsigned int N);

   //   FORWARD FOURIER TRANSFORM, INPLACE VERSION
   //     Data - both input data and output
   //     N    - length of both input data and result
   static bool Forward(complex *const Data, const unsigned int N);

   //   INVERSE FOURIER TRANSFORM
   //     Input  - input data
   //     Output - transform result
   //     N      - length of both input data and result
   //     Scale  - if to scale result
   static bool Inverse(const complex *const Input, complex *const Output,
      const unsigned int N, const bool Scale = true);

   //   INVERSE FOURIER TRANSFORM, INPLACE VERSION
   //     Data  - both input data and output
   //     N     - length of both input data and result
   //     Scale - if to scale result
   static bool Inverse(complex *const Data, const unsigned int N,
      const bool Scale = true);

protected:
   //   Rearrange function and its inplace version
   static void Rearrange(const complex *const Input, complex *const Output,
      const unsigned int N);
   static void Rearrange(complex *const Data, const unsigned int N);

   //   FFT implementation
   static void Perform(complex *const Data, const unsigned int N,
      const bool Inverse = false);

   //   Scaling of inverse FFT result
   static void Scale(complex *const Data, const unsigned int N);
};

Class has four public methods for performing FFT: two functions for forward transform and two ones for inverse transform. Every couple consists of in-place version and version that preserves input data and outputs transform result into provided array.

Protected section of the class has as well four functions: two functions for data preprocessing — putting them into convenient order, core function for transform performing and auxiliary function for scaling the result of inverse transform.

Every of four public methods is very similar and is indeed wrapper that controls processing workflow. Let us see how one of them designed.

//   FORWARD FOURIER TRANSFORM, INPLACE VERSION
//     Data - both input data and output
//     N    - length of both input data and result
bool CFFT::Forward(complex *const Data, const unsigned int N)
{
   //   Check input parameters
   if (!Data || N < 1 || N & (N - 1))
      return false;
   //   Rearrange
   Rearrange(Data, N);
   //   Call FFT implementation
   Perform(Data, N);
   //   Succeeded
   return true;
}

Inside wrapper you can find check of input parameters, then data preparation — rearrangement, — and transform itself. Rearrange function uses our “mirrored mathematics” to define new position for every element and swaps elements:

//   Inplace version of rearrange function
void CFFT::Rearrange(complex *const Data, const unsigned int N)
{
   //   Swap position
   unsigned int Target = 0;
   //   Process all positions of input signal
   for (unsigned int Position = 0; Position < N; ++Position)
   {
      //   Only for not yet swapped entries
      if (Target > Position)
      {
         //   Swap entries
         const complex Temp(Data[Target]);
         Data[Target] = Data[Position];
         Data[Position] = Temp;
      }
      //   Bit mask
      unsigned int Mask = N;
      //   While bit is set
      while (Target & (Mask >>= 1))
         //   Drop bit
         Target &= ~Mask;
      //   The current bit is 0 - set it
      Target |= Mask;
   }
}

While loop implements addition of 1 in mirrored manner to get target position for the element.

Now there is turn of the core method, which performs our fast Fourier transform.

//   FFT implementation
void CFFT::Perform(complex *const Data, const unsigned int N,
   const bool Inverse /* = false */)
{
   const double pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846;
   //   Iteration through dyads, quadruples, octads and so on...
   for (unsigned int Step = 1; Step < N; Step <<= 1)
   {
      //   Jump to the next entry of the same transform factor
      const unsigned int Jump = Step << 1;
      //   Angle increment
      const double delta = pi / double(Step);
      //   Auxiliary sin(delta / 2)
      const double Sine = sin(delta * .5);
      //   Multiplier for trigonometric recurrence
      const complex Multiplier(-2. * Sine * Sine, sin(delta));
      //   Start value for transform factor, fi = 0
      complex Factor(1.);
      //   Iteration through groups of different transform factor
      for (unsigned int Group = 0; Group < Step; ++Group)
      {
         //   Iteration within group 
         for (unsigned int Pair = Group; Pair < N; Pair += Jump)
         {
            //   Match position
            const unsigned int Match = Pair + Step;
            //   Second term of two-point transform
            const complex Product(Factor * Data[Match]);
            //   Transform for fi + pi
            Data[Match] = Data[Pair] - Product;
            //   Transform for fi
            Data[Pair] += Product;
         }
         //   Successive transform factor via trigonometric recurrence
         Factor = Multiplier * Factor + Factor;
      }
   }
}

The code is exact reflection of our FFT algorithm and butterfly scheme in fig. 2. Difference between forward and inverse transforms is reflected in the first line of the method where the proper sign for exponent argument is set. Initializations inside outer loop are just preparations for successive calculation of factors via trigonometric recurrence. And the job done inside inner loop, which performs butterfly operations. Trigonometric recurrence in the last line is exactly our expression (10).

Wrappers for inverse transform are designed the same way as for forward one:

//   INVERSE FOURIER TRANSFORM, INPLACE VERSION
//     Data  - both input data and output
//     N     - length of both input data and result
//     Scale - if to scale result
bool CFFT::Inverse(complex *const Data, const unsigned int N,
   const bool Scale /* = true */)
{
   //   Check input parameters
   if (!Data || N < 1 || N & (N - 1))
      return false;
   //   Rearrange
   Rearrange(Data, N);
   //   Call FFT implementation
   Perform(Data, N, true);
   //   Scale if necessary
   if (Scale)
      CFFT::Scale(Data, N);
   //   Succeeded
   return true;
}

The only difference is conditional scaling of the result at postprocessing stage. By default the scaling is performed according to (13) but if one interested only in relative values she can drop the corresponding flag to skip this operation. Scaling itself is a primitive function below.

//   Scaling of inverse FFT result
void CFFT::Scale(complex *const Data, const unsigned int N)
{
   const double Factor = 1. / double(N);
   //   Scale all data entries
   for (unsigned int Position = 0; Position < N; ++Position)
      Data[Position] *= Factor;
}

Well, our FFT is ready. You can download full source code here:

Download fast Fourier transform — FFT — C++ source code (zip, 4 Kb)

Full file listings are available online as well:

Optimized for high performance source code of complex number class you can find here:

6. How to use

To utilize the FFT you should include header file fft.h and place in your code lines like:

...Usually at the top of the file
//   Include FFT header
#include "fft.h"
.
...Some code here
.
   ...Inside your signal processing function
   //   Allocate memory for signal data
   complex *pSignal = new complex[1024];
   ...Fill signal array with data
   //   Apply FFT
   CFFT::Forward(pSignal, 1024);
   ...Utilize transform result
   //   Free memory
   delete[] pSignal;

Here in-place forward Fourier transform performed for signal of 1024 samples size.

7. Final remarks

The FFT algorithm implemented in literature is called Cooley-Tukey. There is also Sand-Tukey algorithm that rearranges data after performing butterflies and in its case butterflies look like ours in fig. 2 but mirrored to the right so that big butterflies come first and small ones do last.

From all our considerations follows that length of input data for our algorithm should be power of 2. In the case length of the input data is not power of 2 it is good idea to extend the data size to the nearest power of 2 and pad additional space with zeroes or input signal itself in a periodic manner — just copy necessary part of the signal from its beginning. Padding with zeroes usually works well.

If you were attentive you could notice that butterflies in parenthesis and brackets in (7) do not really need multiplications but only additions and subtractions. So, optimizing two deepest levels of butterflies we can even improve FFT performance.

Posted in Image / Video Filters | Tagged: , | Leave a Comment »

Gaussian filter, or Gaussian blur

Posted by Hemprasad Y. Badgujar on October 12, 2014


Gaussian filter, or Gaussian blur

Category. Digital signal and image processing (DSP and DIP) software development.

Abstract. The article is a practical tutorial for Gaussian filter, or Gaussian blur understanding and implementation of its separable version. Article contains theory, C++ source code, programming instructions and sample application.

Download Gaussian filter, or Gaussian blur for Win32 (zip, 890 Kb)

Download Gaussian filter, or Gaussian blur C++ source code (zip, 4 Kb)

Original and blurred images

1. Introduction to Gaussian filter, or Gaussian blur

Gaussian filter is windowed filter of linear class, by its nature is weighted mean. Named after famous scientist Carl Gauss because weights in the filter calculated according to Gaussian distribution — the function Carl used in his works. Another name for this filter is Gaussian blur.

To get acquainted with filter window idea in signal and image processing read our “Filter window, or filter mask” article.

2. Understanding Gaussian filter, or Gaussian blur

First of all let us have a look at what that Gaussian distribution is. Gaussian distribution, or normal distribution, is really a function of probability theory. Often this function is referenced as bell-function because of its shape. The most general function formula is:

Gaussian distribution (1)

And its plot is depicted below — fig. 1.

Fig. 1. Gaussian or normal distribution.Fig. 1. Gaussian or normal distribution.

In our case we can suppose parameter a — which called distribution mean or statistical expectation — responsible for distribution shifting along x axis to be zero: a=0; and work with simplified form:

Gaussian distribution (2)

Thus the function is negative exponential one of squared argument. Argument divider σ plays the role of scale factor. σ parameter has special name: standard deviation, and its square σ2variance. Premultiplier in front of the exponent is selected the area below plot to be 1. Pay attention the function is defined everywhere on real axis x∈(-∞, ∞) which means it spreads endlessly to the left and to the right.

Now, first point is we are working in discrete realm, so our Gaussian distribution turns into the set of values at discrete points.

Second, we cannot work with something that spreads endlessly to the left and to the right. It means Gaussian distribution is to be truncated. The question is — where? Usually in practice used the rule of 3σ that is the part of Gaussian distribution utilized is x∈[-3σ, 3σ] — see fig. 2.

Fig. 2. Gaussian distribution truncated at points ±3σ.Fig. 2. Gaussian distribution truncated at points ±3σ.

Why? Good question. To understand that let us see how much we have trimmed. The area below truncated part is:

Area below truncated Gaussian distribution (3)

Now we can start up MatLab and type lines like:

G=inline('1/sqrt(2*pi)*exp(-x.^2/2)');
quad(G, -3, 3)

Which tells MatLab to calculate our integral (3). The result is

ans =
    0.9973

So, the area below our trimmed part is ≈0.997 — that is function is trimmed at points where we have accuracy better than 0.5%. Very good in our case.

Back to our discrete normal distribution: we are interested in points {x-N, xN+1, … , x-1, x0, x1, … , xN-1, xN}, where x-n=-xn and xN=3σ, and respectively in values of normal distribution at these points: {G(x-N), G(xN+1), … , G(x-1), G(x0), G(x1), … , G(xN-1), G(xN)}. So, we have 2N+1 value set {Gn | n=-N, –N+1, … , N} where Gn=G(xn).

What shall we do now with this set of values? Use it as window weights! That means we have weighted window of 2N+1 size. To be perfect we are to scale our Gn: G’n=Gn/k so that ∑G’n=1 — sum of all of them to be one: G’n=Gn/∑Gn — which means k=∑Gn.

Thus, we have our window weights {G’n | n=-N, –N+1, … , N} where G’n=Gn/k, Gn=G(xn), k=∑Gn and xN=3σ which means as well xn=3σn/N. To get expression for G’n for practical use let us write down the detailed formula:

Weights calculation (4)

As you can see we have simplification: σ is eliminated from consideration and G’n could be calculated easier via new values G”n and their sumk’.

How to use these weights? If we have some input signal S={si} then for every signal element si new modified value s’i will be s’i=∑G’nsi+n, n=-N, –N+1, … , N. In words that means “for every element put our window so that this element is in the center of the window, multiply every element in the window by corresponding weight and sum up all those products, the sum got is new filtered value”.

Now we can write down step-by-step instructions for processing by 1D Gaussian filter or blur.

1D Gaussian filter, or Gaussian blur algorithm:

  1. Given window size 2N+1 calculate support points xn=3n/N, n=-N, –N+1, … , N;
  2. Calculate values G”n;
  3. Calculate scale factor k’=∑G”n;
  4. Calculate window weights G’n=G”n/k’;
  5. For every signal element:
    1. Place window over it;
    2. Pick up elements;
    3. Multiply elements by corresponding window weights;
    4. Sum up products — this sum is new filtered value.

Let us proceed with 2D case.

3. 2D case

Expression for 2D Gaussian distribution is:

2D Gaussian distribution (5)

And its plot is depicted below — fig. 3.

Fig. 3. 2D Gaussian or normal distribution.Fig. 3. 2D Gaussian or normal distribution.

Gaussian distribution has surprising property. Look, its expression could be rewritten as:

2D Gaussian distribution separation (6)

Which means 2D distribution is split into a pair of 1D ones, and so, 2D filter window (fig. 3) is separated into a couple of 1D ones from fig. 2. Filter version that utilizes this fact is called separable one.

In practice it means that to apply filter to an image it is enough to filter it in horizontal direction with 1D filter and then filter the result with the same filter in vertical direction. Which direction first really makes no difference — our operation is commutative. Thus,

2D separable Gaussian filter, or Gaussian blur, algorithm:

  1. Calculate 1D window weights G’n;
  2. Filter every image line as 1D signal;
  3. Filter every filtered image column as 1D signal.

2D Gaussian filtering with [2N+1]×[2N+1] window is reduced to a couple of 1D filterings with 2N+1 window. That means significant speed-up especially for large images because of jump from O(N2) to O(N) number of operations.

Now, when we have the algorithm, it is time to write some code — let us come down to programming.

4. 2D Gaussian filter, or 2D Gaussian blur programming

We are starting with 2D filter because 1D one could be easily got just by treating signal as one-line image and canceling vertical filtering.

First of all a couple of simple auxiliary structures.

//   Internal auxiliary structure - data array size descriptor
struct CSize
{
   unsigned int x;   //   Array width
   unsigned int y;   //   Array height

   //   Default constructor
   CSize(): x(0), y(0) {}
   //   Constructor with initialization
   CSize(unsigned int _x, unsigned int _y): x(_x), y(_y) {}

   //   Initialization
   void Set(unsigned int _x, unsigned int _y) { x = _x; y = _y; }
   //   Area
   unsigned int Area() const { return x * y; }
};

//   Internal auxiliary structure - array descriptor
struct CArray
{
   CSize Size;   //   Array size
   T *Buffer;    //   Element buffer

   //   Default constructor
   CArray(): Buffer(NULL) {}
   //   Constructors with initialization
   CArray(T *_Buffer, const CSize &_Size): Buffer(_Buffer), Size(_Size) {}
   CArray(T *_Buffer, unsigned int _N): Buffer(_Buffer), Size(_N, 1) {}
};

CSize structure keeps the image size and CArray structure is image descriptor — its size and pointer to image data. Now more complicated structure — our Gaussian window.

//   Internal auxiliary structure - filter window descriptor
struct CWindow
{
   double *Weights;   //   Window weights
   unsigned int Size;   //   Window size

   //   Default constructor
   CWindow(): Weights(NULL), Size(0), Correction(.5 - double(T(.5))) {}
   //   Destructor
   ~CWindow() { if (Weights) delete[] Weights; }

   //   Filter window creation
   bool Create(unsigned int _Size);

   //   FILTER WINDOW APPLICATION
   //     _Element - start element in signal/image
   T Apply(const T *_Element) const
   {
      //   Apply filter - calculate weighted mean
      double Sum = 0.;
      const double *WeightIter = Weights;
      const T *ElIter = _Element;
      const double *const End = Weights + Size;
      while (WeightIter < End)
         Sum += *(WeightIter++) * double(*(ElIter++));
      return T(Sum + Correction);
   }

protected:
   const double Correction;   //   Result correction
};

CWindow structure designed to keep window size and set of weights, as well it has two methods to calculate weights and to apply 1D Gaussian filter starting from a given image element: Create and Apply. As you can see Apply code is trivial and code for Create is a little bit more complicated:

//   FILTER WINDOW CREATION
//     _Size - window size
template <class T> bool TGaussianBlur<T>::CWindow::Create(unsigned int _Size)
{
   //   Allocate window buffer
   Weights = new double[_Size];
   //   Check allocation
   if (!Weights)
      //   Window buffer allocation failed
      return false;
   //   Set window size
   Size = _Size;
   //   Window half
   const unsigned int Half = Size >> 1;
   //   Central weight
   Weights[Half] = 1.;
   //   The rest of weights
   for (unsigned int Weight = 1; Weight < Half + 1; ++Weight)
   {
      //   Support point
      const double x = 3.* double(Weight) / double(Half);
      //   Corresponding symmetric weights
      Weights[Half - Weight] = Weights[Half + Weight] = exp(-x * x / 2.);
   }
   //   Weight sum
   double k = 0.;
   for (unsigned int Weight = 0; Weight < Size; ++Weight)
      k += Weights[Weight];
   //   Weight scaling
   for (unsigned int Weight = 0; Weight < Size; ++Weight)
      Weights[Weight] /= k;
   //   Succeeded
   return true;
}

The method implements steps 1–4 of 1D Gaussian filter algorithm to calculate weights according to expression (4) and utilizes window symmetry G-n=Gn.

Now, the last problem to be solved before we can start filtering the image is its extension.

5. Extension

There is a problem every windowed filter deals with. Being placed at the edge filter window lacks for data elements to be processed. There are two ways to solve the problem: first is not to process edges and second, cleverer one, to extend data across edges. We are taking the second approach and extending our data like depicted in fig. 4.

Fig. 4. Data extension.Fig. 4. Data extension.

And here is CExtension class that will do that job for us.

//   Internal auxiliary structure - array extension descriptor
struct CExtension: public CArray
{
   unsigned int Margin;   //   Extension margins

   enum EMode {ModeHorizontal, ModeVertical};

   //   Default cosntructor
   CExtension(): Margin(0), Mode(ModeHorizontal) {}
   //   Destructor
   ~CExtension() { if (Buffer) delete[] Buffer; }

   //   Mode setting
   void SetMode(EMode _Mode) { Mode = _Mode; }
   //   Extension memory allocation
   bool Allocate(unsigned int _N, unsigned int _W)
   { return _Allocate(CSize(_N, 1), _W >> 1); }
   bool Allocate(const CSize &_Size, unsigned int _W)
   { return _Allocate(_Size, _W >> 1); }
   //   Pasting data into extension from data array
   void Paste(const T * _Start);
   //   Extension
   void Extend();

protected:
   EMode Mode;   //   Current mode

   //   Extension memory allocation
   bool _Allocate(const CSize &_Size, unsigned int _Margin);
};

The job is done inside methods _Allocate, Paste and Extend. _Allocate method allocates space in memory enough to keep extended image line or column.

//   EXTENSION MEMORY ALLOCATION
//     _Size   - signal/image size
//     _Margin - extension margins
template <class T> bool TGaussianBlur<T>::CExtension::_Allocate(
   const CSize &_Size, unsigned int _Margin)
{
   //   Allocate extension buffer
   Buffer = new T[(_Size.x > _Size.y ? _Size.x : _Size.y) + (_Margin << 1)];
   //   Check buffer allocation
   if (!Buffer)
      //   Buffer allocation failed
      return false;
   //   Initialize size descriptors
   Size = _Size;
   Margin = _Margin;
   //   Succeeded
   return true;
}

Method Paste inserts image line or column at proper place in extension.

//   PASTING DATA LINE/COLUMN INTO EXTENSION FROM DATA ARRAY
//     _Start - start postion in image/signal to paste from
template <class T> void TGaussianBlur<T>::CExtension::Paste(const T *const _Start)
{
   if (Mode == ModeHorizontal)
   {
      //   Paste line
      memcpy(Buffer + Margin, _Start, Size.x * sizeof(T));
   }
   else
   {
      //   Stop position
      const T *const Stop = _Start + Size.Area();
      //   Array iterator
      const T *ArrIter = _Start;
      //   Extension iterator
      T *ExtIter = Buffer + Margin;
      //   Paste array column element by element
      while (ArrIter < Stop)
      {
         //   Copy line
         *(ExtIter++) = *ArrIter;
         //   Jump to the next line
         ArrIter += Size.x;
      }
   }
}

Pasting has two modes — horizontal and vertical. In horizontal mode line is copied into the extension and in vertical mode column is copied element by element. Method Extend mirrors pasted data.

//   EXTENSION CREATION
template <class T> void TGaussianBlur<T>::CExtension::Extend()
{
   //   Line size
   const unsigned int Line = Mode == ModeHorizontal ? Size.x : Size.y;
   //   Stop position
   const T *const Stop = Buffer - 1;
   //   Left extension iterator
   T *ExtLeft = Buffer + Margin - 1;
   //   Left array iterator
   const T *ArrLeft = ExtLeft + 2;
   //   Right extension iterator
   T *ExtRight = ExtLeft + Line + 1;
   //   Left array iterator
   const T *ArrRight = ExtRight - 2;
   //   Create extension line element by element
   while (ExtLeft > Stop)
   {
      //   Left extension
      *(ExtLeft--) = *(ArrLeft++);
      //   Right extension
      *(ExtRight++) = *(ArrRight--);
   }
}

Now we have everything to program filtering method.

//   2D GAUSSIAN BLUR
//     pImage  - input image
//     pResult - output image, NULL for inplace processing
//     N       - width of the image
//     M       - height of the image
//     W       - window size
template <class T> bool TGaussianBlur<T>::Filter(T *pImage, T *pResult,
   unsigned int N, unsigned int M, unsigned int W) const
{
   //   Check input data consistency
   if (!Consistent(pImage, CSize(N, M), W))
     return false;
   //   Allocate extension
   CExtension Extension;
   if (!Extension.Allocate(CSize(N, M), W))
      return false;
   //   Create image descriptor
   CArray Image(pImage, CSize(N, M));
   //   Create filter window
   CWindow Window;
   if (!Window.Create(W))
      return false;
   //   Stop postion
   const T * ExtStop = Extension.Buffer + Extension.Size.x;
   //   Result iterator
   T *ResIter = pResult ? pResult : pImage;
   //   Image iterator
   const T *ImIter = Image.Buffer;
   //   Image stop position
   const T * ImStop = Image.Buffer + Image.Size.Area();
   //   Filter line by line
   while (ImIter < ImStop)
   {
      //   Paste image line into extension
      Extension.Paste(ImIter);
      //   Extend image line
      Extension.Extend();
      //   Extension iterator
      const T *ExtIter = Extension.Buffer;
      //   Apply filter to every pixel of the line
      while (ExtIter < ExtStop)
         *(ResIter++) = Window.Apply(ExtIter++);
         //   Move to the next line
      ImIter += Image.Size.x;
   }
   //   Initialize image descriptor with filter result
   Image.Buffer = pResult ? pResult : pImage;
   //   Set vertical extension mode
   Extension.SetMode(CExtension::ModeVertical);
   //   Extension stop position
   ExtStop = Extension.Buffer + Extension.Size.y;
   //   Result column iterator
   T *ResColumnIter = pResult ? pResult : pImage;
   //   Image iterator
   ImIter = Image.Buffer;
   //   Image stop position
   ImStop = Image.Buffer + Image.Size.x;
   //   Filter column by column
   while (ImIter < ImStop)
   {
      //   Paste image column into extension
      Extension.Paste(ImIter++);
      //   Extend image column
      Extension.Extend();
      //   Extension iterator
      const T *ExtIter = Extension.Buffer;
      //   Result pixel iterator
      ResIter = ResColumnIter;
      //   Apply fitler to every pixel of the column
      while (ExtIter < ExtStop)
      {
         *ResIter = Window.Apply(ExtIter++);
         ResIter += Image.Size.x;
      }
      //   Move to the next column
      ++ResColumnIter;
   }
   //   Succeeded
   return true;
}

Structure of the method is quite straightforward: input parameters check, memory allocation, window weights calculation, applying 1D filter in horizontal direction (first loop) and applying it in vertical direction (second loop). Parameters check is a short function below.

//   Internal auxiliary functions - check input data consistency
bool Consistent(const T *_Image, const CSize &_Size, unsigned int _W) const
{
   return  _Image && _Size.x && _Size.y && _W &&
      _Size.x > (_W >> 1) && _Size.y > (_W >> 1) && _W & 1;
}

Which means pointer to image should not be NULL, image and filter window sizes should be some positive numbers, window size 2N+1≡_Wshould be odd number and N should be less than image size in any direction.

6. 1D Gaussian filter, or Gaussian blur programming

1D filter is truncated version of 2D filter:

//   1D GAUSSIAN BLUR
//     pSignal - input signal;
//     pResult - output signal, NULL for inplace processing
//     N       - length of the signal
//     W       - window size, odd positive number
template <class T> bool TGaussianBlur<T>::Filter(T *pSignal, T *pResult,
   unsigned int N, unsigned int W) const
{
   //   Check input data cosnsitency
   if (!Consistent(pSignal, N, W))
      return false;
   //   Allocate extension
   CExtension Extension;
   if (!Extension.Allocate(N, W))
      return false;
   //   Create signal descriptor
   const CArray Signal(pSignal, N);
   //   Paste signal into extension
   Extension.Paste(Signal.Buffer);
   //   Extend signal
   Extension.Extend();
   //   Create filter window
   CWindow Window;
   if (!Window.Create(W))
      return false;
   //   Extension iterator
   const T *ExtIter = Extension.Buffer;
   //   Extension stop position
   const T *const ExtStop = Extension.Buffer + Extension.Size.x;
   //   Result iterator
   T *ResIter = pResult ? pResult : pSignal;
   //   Filter - apply to every element
   while (ExtIter < ExtStop)
      *(ResIter++) = Window.Apply(ExtIter++);
   //   Succeeded
   return true;
}

You can see familiar steps here just now filter makes one pass along signal array.

Our separable Gaussian filter library is ready! You can download full source code here:

Download Gaussian filter, or Gaussian blur C++ source code (zip, 4 Kb)

Full file listings are available online as well:

6. How to use

Pay attention that our solution is universal. First, our filter has window of arbitrary size. Second, we have developed class template that suits any type of input data —

TGaussianBlur<unsigned char> BlurFilter;
TGaussianBlur<short> BlurFilter;
TGaussianBlur<float> BlurFilter;
TGaussianBlur<int> BlurFilter;
TGaussianBlur<unsigned int> BlurFilter;

— all these declarations are valid and will process your data graciously.

To use the filter you should include header file gaussianblur.h and place in your code lines like:

...Usually at the top of the file
//   Include Gaussian blur header
#include "gaussianblur.h"
.
...Some code here
.
   ...Inside your signal/image processing function
   //   Create instance of Gaussian blur filter
   TGaussianBlur<double> BlurFilter;
   //   Apply Gaussian blur
   BlurFilter.Filter(pImage, NULL, 512, 512, 13);

Here image is stored as double array of 512×512 size and it is filtered in-place with Gaussian window of 13×13 pixels width.

And now — an application to play around!

7. Color Gaussian blur

Download Gaussian filter, or Gaussian blur for Win32 (zip, 890 Kb)

We have created an application to see Gaussian filter in action. The sample package includes 3 files — the applications, sample image and description:

  • gaussianblur.exe — Gaussian filter,
  • sample.bmp — 24-bit sample image,
  • readme.txt — description.

Be aware of the fact, that this sample uses OpenGL, so it should be supported by your system (usually that is the case).

8. How to use

Start up gaussianblur.exe application. Load the image.

Fig. 5. Original image. Screenshot.Fig. 5. Original image.

Set filter window size: Set >> Window size… or click w button in toolbar, in dialog key in size, for instance 13. Blur image with Gaussian filter by choosing Set >> Filter or clicking F-button in toolbar. See the result.

Fig. 6. Image blurred by Gaussian filter. Screenshot.Fig. 6. Image blurred by Gaussian filter.

Posted in Image / Video Filters | Leave a Comment »

Hybrid median filter

Posted by Hemprasad Y. Badgujar on October 12, 2014


Hybrid median filter

Category. Digital image processing (DIP) software development.

Abstract. The article is a practical tutorial for hybrid median filter understanding and implementation. Article contains theory, C++ source code, programming instructions and sample applications.

Download hybrid median filter for Win32 (zip, 1.03 Mb)

Download hybrid median filter C++ source code (zip, 2 Kb)

Corrupted and restored images

1. Introduction to hybrid median filter

Hybrid median filter is windowed filter of nonlinear class, that easily removes impulse noise while preserving edges. In comparison with basic version of the median filter hybrid one has better corner preserving characteristics. The basic idea behind filter is for any element of the signal (image) apply median technique several times varying window shape and then take the median of the got median values.

To get acquainted with filter window idea in signal and image processing read our “Filter window, or filter mask” article.

2. Understanding hybrid median filter

Now let us see, how hybrid version of the median filter works. In one phrase idea is like “apply median filter with cross-mask, apply median filter with x-mask and take the median of got results and element itself”. Hybrid median filter workflow depicted below in fig. 1.

Fig. 1. Hybrid median filter workflow.Fig. 1. Hybrid median filter workflow.

Due to the nature of the workflow the filter dimension is not less than 2D.

Now let us write down step-by-step instructions for processing by hybrid median filter.

Hybrid median filter algorithm:

  1. Place a cross-window over element;
  2. Pick up elements;
  3. Order elements;
  4. Take the middle element;
  5. Place a x-window over element;
  6. Pick up elements;
  7. Order elements;
  8. Take the middle element;
  9. Pick up result in point 4, 8 and element itself;
  10. Order elements;
  11. Take the middle element.

Now, when we have the algorithm, it is time to write some code — let us come down to programming.

3. Hybrid median filter programming

We can see that operation of taking median — “order elements” and “take the middle element” — repeated three times in the list. So it is good idea to put those steps into separate function:

//   MEDIAN calculation
//     elements - input elements
//     N        - number of input elements
element median(element* elements, int N)
{
   //   Order elements (only half of them)
   for (int i = 0; i < (N >> 1) + 1; ++i)
   {
      //   Find position of minimum element
      int min = i;
      for (int j = i + 1; j < N; ++j)
         if (elements[j] < elements[min])
            min = j;
      //   Put found minimum element in its place
      const element temp = elements[i];
      elements[i] = elements[min];
      elements[min] = temp;
   }
   //   Get result - the middle element
   return elements[N >> 1];
}

In the function above we are using a code optimization trick. Since everything we need is middle element, it is enough to order just a half of elements.

Let us have 2D signal or image of length N×M as input and let we choose filter window of 3×3 size. The first step is window placing — we do that by changing index of the leading element:

//   Move window through all elements of the image
for (int m = 1; m < M - 1; ++m)
   for (int n = 1; n < N - 1; ++n)

Pay attention, that we are not processing first and last rows and columns. The problem is we cannot process elements at edges, because in this case the part of the filter window is empty. We will discuss below, how to solve that problem.

The second step is picking up elements around with cross-window, ok:

//   Pick up cross-window elements
window[0] = image[(m - 1) * N + n];
window[1] = image[m * N + n - 1];
window[2] = image[m * N + n];
window[3] = image[m * N + n + 1];
window[4] = image[(m + 1) * N + n];

The third and fourth steps are done by our auxiliary function:

//   Get median
results[0] = median(window, 5);

Fifth and sixth steps are picking up elements with x-window:

//   Pick up x-window elements
window[0] = image[(m - 1) * N + n - 1];
window[1] = image[(m - 1) * N + n + 1];
window[2] = image[m * N + n];
window[3] = image[(m + 1) * N + n - 1];
window[4] = image[(m + 1) * N + n + 1];

Seventh and eighth steps are made again by our auxiliary function:

//   Get median
results[1] = median(window, 5);

To complete ninth step enough to pick up our leading element itself:

//   Pick up leading element
results[2] = image[m * N + n];

Final steps are performed by our useful function:

//   Get result
result[(m - 1) * (N - 2) + n - 1] = median(results, 3);

At last, let us write down the entire algorithm as function:

//   2D HYBRID MEDIAN FILTER implementation
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
void _hybridmedianfilter(const element* image, element* result, int N, int M)
{
   //   Move window through all elements of the image
   for (int m = 1; m < M - 1; ++m)
      for (int n = 1; n < N - 1; ++n)
      {
         element window[5];
         element results[3];
         //   Pick up cross-window elements
         window[0] = image[(m - 1) * N + n];
         window[1] = image[m * N + n - 1];
         window[2] = image[m * N + n];
         window[3] = image[m * N + n + 1];
         window[4] = image[(m + 1) * N + n];
         //   Get median
         results[0] = median(window, 5);
         //   Pick up x-window elements
         window[0] = image[(m - 1) * N + n - 1];
         window[1] = image[(m - 1) * N + n + 1];
         window[2] = image[m * N + n];
         window[3] = image[(m + 1) * N + n - 1];
         window[4] = image[(m + 1) * N + n + 1];
         //   Get median
         results[1] = median(window, 5);
         //   Pick up leading element
         results[2] = image[m * N + n];
         //   Get result
         result[(m - 1) * (N - 2) + n - 1] = median(results, 3);
      }
}

Looks very straightforward.

Type element could be defined as:

typedef double element;

4. Treating edges

For all window filters there is some problem. That is edge treating. If you place window over an element at the edge, some part of the window will be empty. To fill the gap, signal should be extended. For hybrid median filter there is good idea to extend image symmetrically, like this:

Fig. 2. Image extension.Fig. 2. Image extension.

In other words we are adding lines at the top and at the bottom of the image and add columns to the left and to the right of it.

So, before passing signal to our median filter function the signal should be extended. Let us write down the wrapper, that makes all preparations.

//   2D HYBRID MEDIAN FILTER wrapper
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
void hybridmedianfilter(element* image, element* result, int N, int M)
{
   //   Check arguments
   if (!image || N < 1 || M < 1)
      return;
   //   Allocate memory for signal extension
   element* extension = new element[(N + 2) * (M + 2)];
   //   Check memory allocation
   if (!extension)
      return;
   //   Create image extension
   for (int i = 0; i < M; ++i)
   {
      memcpy(extension + (N + 2) * (i + 1) + 1,
         image + N * i,
         N * sizeof(element));
      extension[(N + 2) * (i + 1)] = image[N * i];
      extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1];
   }
   //   Fill first line of image extension
   memcpy(extension,
      extension + N + 2,
      (N + 2) * sizeof(element));
   //   Fill last line of image extension
   memcpy(extension + (N + 2) * (M + 1),
      extension + (N + 2) * M,
      (N + 2) * sizeof(element));
   //   Call hybrid median filter implementation
   _hybridmedianfilter(extension, result ? result : image, N + 2, M + 2);
   //   Free memory
   delete[] extension;
}

As you can see, our code takes into account some practical issues. First of all we check our input parameters — image should not be NULL, and image sizes should be positive:

//   Check arguments
if (!image || N < 1 || M < 1)
   return;

Now let us allocate memory for signal extension.

//   Allocate memory for signal extension
element* extension = new element[(N + 2) * (M + 2)];

And check memory allocation.

//   Check memory allocation
if (!extension)
   return;

Now we are building extension.

//   Create signal extension
for (int i = 0; i < M; ++i)
{
   memcpy(extension + (N + 2) * (i + 1) + 1,
      image + N * i,
      N * sizeof(element));
   extension[(N + 2) * (i + 1)] = image[N * i];
   extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1];
}
//   Fill first line of image extension
memcpy(extension,
   extension + N + 2,
   (N + 2) * sizeof(element));
//   Fill last line of image extension
memcpy(extension + (N + 2) * (M + 1),
   extension + (N + 2) * M,
   (N + 2) * sizeof(element));

Finally, everything is ready — filtering!

//   Call hybrid median filter implementation
_hybridmedianfilter(extension, result ? result : image, N + 2, M + 2);

And to complete the job — free memory.

//   Free memory
delete[] extension;

Since we are using memory management function from standard library, we should include its header.

#include <memory.h>

5. Hybrid median filter declaration file

Now we have three functions: helper function, hybrid median filter and entry point that makes preparations. Let us write code for header file.

#ifndef _HYBRIDMEDIANFILTER_H_
#define _HYBRIDMEDIANFILTER_H_

//   Image element type
typedef double element;

//   2D HYBRID MEDIAN FILTER, window size 3x3
//     image  - input image
//     result - output image, NULL for inplace processing
//     N      - width of the image
//     M      - height of the image
void hybridmedianfilter(element* image, element* result, int N, int M);

#endif

Now all is ready. The code we have written is good both for Linux/Unix and Windows platforms. You can download full hybrid median filter source code here:

Download hybrid median filter C++ source code (zip, 2 Kb)

Full file listings are available online as well:

And now — a couple of applications to play around!

6. Color median filter: image restoration

Download hybrid median filter for Win32 (zip, 1.03 Mb)

We have created a couple of applications to show hybrid median filter capabilities in restoration images corrupted by impulse noise. The sample package includes 4 files — two applications, sample image and description:

  • hybridmedian.exe — median filter,
  • corrupter.exe — destructive noise generator,
  • sample.bmp — 24-bit sample image,
  • readme.txt — description.

Be aware of the fact, that this sample uses OpenGL, so it should be supported by your system (usually that is the case).

7. Step 1: prepare corrupted image

We have created impulse noise generator that will help us to prepare corrupted image. Start up corrupter.exe application and load image to be corrupted. Choose Set >> Corruption… or click N button in toolbar and set noise level at 5–15%. Click OK. Then save corrupted image.

Fig. 3. Corruption by impulse noise. Screenshot.Fig. 3. Corruption by impulse noise.

8. Step 2: restore corrupted image

Start up hybridmedian.exe application. Load the saved corrupted image. Apply hybrid median filter by choosing Set >> Filter or clickingF-button in toolbar. See the result. If necessary, filter the image one more time.

Fig. 4. Image restored by hybrid median filter. Screenshot.Fig. 4. Image restored by hybrid median filter.

 

Posted in Image / Video Filters | Tagged: | Leave a Comment »

Alpha-trimmed mean filter

Posted by Hemprasad Y. Badgujar on October 12, 2014


Alpha-trimmed mean filter

Category. Digital signal and image processing (DSP and DIP) software development.

Abstract. The article is a practical tutorial for alpha-trimmed mean filter understanding and implementation. Article contains theory, C++ source code, programming instructions and sample applications.

Download alpha-trimmed mean filter for Win32 (zip, 603 Kb)

Download alpha-trimmed mean filter C++ source code (zip, 2 Kb)

Corrupted and restored images

1. Introduction to alpha-trimmed mean filter

Alpha-trimmed mean filter is windowed filter of nonlinear class, by its nature is hybrid of the mean and median filters. The basic idea behind filter is for any element of the signal (image) look at its neighborhood, discard the most atypical elements and calculate mean value using the rest of them. Alpha you can see in the name of the filter is indeed parameter responsible for the number of trimmed elements.

To get acquainted with filter window idea in signal and image processing read our “Filter window, or filter mask” article.

2. Understanding alpha-trimmed mean filter

Now let us see, how to get alpha-trimmed mean value in practice. The basic idea here is to order elements, discard elements at the beginning and at the end of the got ordered set and then calculate average value using the rest. For instance, let us calculate alpha-trimmed mean for the case, depicted in fig. 1.

Fig. 1. Alpha-trimmed mean calculation.Fig. 1. Alpha-trimmed mean calculation.

Thus, to get an alpha-trimmed mean we are to order elements, eliminate elements at the beginning and at the end of the got sequenced collection and get average of the remaining elements. That is all — now alpha-trimmed mean calculated and signal, 1D in our case, filtered by alpha-trimmed mean filter! Let us make resume and write down step-by-step instructions for processing by alpha-trimmed mean filter.

Alpha-trimmed mean filter algorithm:

  1. Place a window over element;
  2. Pick up elements;
  3. Order elements;
  4. Discard elements at the beginning and at the end of the got ordered set;
  5. Take an average — sum up the remaining elements and divide the sum by their number.

A couple of words about alpha parameter responsible for trimmed elements. In practice alpha is the number of elements to be discarded, for instance, in our case alpha is two. Since our filter is symmetric one alpha is an even nonnegative number less than size of the filter window. Minimum value for alpha parameter is zero and in this case alpha-trimmed mean filter degenerates into mean filter. Maximum value for alpha is filter window size minus one and in this case filter degenerates into median filter.

Now, when we have the algorithm, it is time to write some code — let us come down to programming.

3. 1D alpha-trimmed mean filter programming

In this section we develop 1D alpha-trimmed mean filter with window of size 5. Let us have 1D signal of length N as input. The first step is window placing — we do that by changing index of the leading element:

//   Move window through all elements of the signal
for (int i = 2; i < N - 2; ++i)

Pay attention, that we are starting with the third element and finishing with the last but two. The problem is we cannot start with the first element, because in this case the left part of the filter window is empty. We will discuss below, how to solve that problem.

The second step is picking up elements around, ok:

//   Pick up window elements
for (int j = 0; j < 5; ++j)
   window[j] = signal[i - 2 + j];

The third step is putting window elements in order. But we will use a code optimization trick here. As soon as we need only middle part of the ordered set we can omit putting in order last α/2 elements. So, it is enough to process just 5-α/2 elements:

//   End of the trimmed ordered set
const int end = 5 - (alpha >> 1);

Now:

//   Order elements (only necessary part of them)
for (int j = 0; j < end; ++j)
{
   //   Find position of minimum element
   int min = j;
   for (int k = j + 1; k < 5; ++k)
      if (window[k] < window[min])
         min = k;
   //   Put found minimum element in its place
   const element temp = window[j];
   window[j] = window[min];
   window[min] = temp;
}

To discard elements at the beginning of the ordered set let us define the proper constant:

//   Start of the trimmed ordered set
const int start = alpha >> 1;

The final step — taking the average:

//   Get result - the mean value of the elements in trimmed set
result[i - 2] = window[start];
for (int j = start + 1; j < end; ++j)
   result[i - 2] += window[j];
result[i - 2] /= N - alpha;

At last, let us write down the entire algorithm as function:

//   1D ALPHA-TRIMMED MEAN FILTER implementation
//     signal - input signal
//     result - output signal
//     N      - length of the signal
//     alpha  - filter alpha parameter
void _alphatrimmedmeanfilter(const element* signal, element* result,
   int N, int alpha)
{
   //   Start of the trimmed ordered set
   const int start = alpha >> 1;
   //   End of the trimmed ordered set
   const int end = 5 - (alpha >> 1);
   //   Move window through all elements of the signal
   for (int i = 2; i < N - 2; ++i)
   {
      //   Pick up window elements
      element window[5];
      for (int j = 0; j < 5; ++j)
         window[j] = signal[i - 2 + j];
      //   Order elements (only necessary part of them)
      for (int j = 0; j < end; ++j)
      {
         //   Find position of minimum element
         int min = j;
         for (int k = j + 1; k < 5; ++k)
            if (window[k] < window[min])
               min = k;
         //   Put found minimum element in its place
         const element temp = window[j];
         window[j] = window[min];
         window[min] = temp;
      }
      //   Get result - the mean value of the elements in trimmed set
      result[i - 2] = window[start];
      for (int j = start + 1; j < end; ++j)
         result[i - 2] += window[j];
      result[i - 2] /= 5 - alpha;
   }
}

Type element could be defined as:

typedef double element;

4. Treating edges

For all window filters there is some problem. That is edge treating. If you place window over first (last) element, the left (right) part of the window will be empty. To fill the gap, signal should be extended. For alpha-trimmed mean filter there is good idea to extend signal symmetrically, like this:

Fig. 2. Signal extension.Fig. 2. Signal extension.

So, before passing signal to our alpha-trimmed mean filter function the signal should be extended. Let us write down the wrapper, that makes all preparations.

//   1D ALPHA-TRIMMED MEAN FILTER wrapper
//     signal - input signal
//     result - output signal
//     N      - length of the signal
//     alpha  - filter alpha parameter
void alphatrimmedmeanfilter(element* signal, element* result, int N, int alpha)
{
   //   Check arguments
   if (!signal || N < 1 || alpha < 0 || 4 < alpha || alpha & 1)
      return;
   //   Treat special case N = 1
   if (N == 1)
   {
      if (result)
         result[0] = signal[0];
      return;
   }
   //   Allocate memory for signal extension
   element* extension = new element[N + 4];
   //   Check memory allocation
   if (!extension)
      return;
   //   Create signal extension
   memcpy(extension + 2, signal, N * sizeof(element));
   for (int i = 0; i < 2; ++i)
   {
      extension[i] = signal[1 - i];
      extension[N + 2 + i] = signal[N - 1 - i];
   }
   //   Call alpha-trimmed mean filter implementation
   _alphatrimmedmeanfilter(extension + 2, result ? result : signal,
      N + 4, alpha);
   //   Free memory
   delete[] extension;
}

As you can see, our code takes into account some practical issues. First of all we check our input parameters — signal should not be NULL, signal length should be positive, alpha parameter should be even nonnegative number less than window size:

//   Check arguments
if (!signal || N < 1 || alpha < 0 || 4 < alpha || alpha & 1)
   return;

Second step — we check case N=1. This case is special one, because to build extension we need at least two elements. For the signal of 1 element length the result is the signal itself. As well pay attention, our alpha-trimmed mean filter works in-place, if output parameter resultis NULL.

//   Treat special case N = 1
if (N == 1)
{
   if (result)
      result[0] = signal[0];
   return;
}

Now let us allocate memory for signal extension.

//   Allocate memory for signal extension
element* extension = new element[N + 4];

And check memory allocation.

//   Check memory allocation
if (!extension)
   return;

Now we are building extension.

//   Create signal extension
memcpy(extension + 2, signal, N * sizeof(element));
for (int i = 0; i < 2; ++i)
{
   extension[i] = signal[1 - i];
   extension[N + 2 + i] = signal[N - 1 - i];
}

Finally, everything is ready — filtering!

//   Call alpha-trimmed mean filter implementation
_alphatrimmedmeanfilter(extension, result ? result : signal, N + 4, alpha);

And to complete the job — free memory.

//   Free memory
delete[] extension;

Since we are using memory management function from standard library, we should include its header.

#include <memory.h>

5. 2D alpha-trimmed mean filter programming

In 2D case we have 2D signal, or image. The idea is the same, just now alpha-trimmed mean filter has 2D window. Window influences only the elements selection. All the rest is the same: ordering elements, trimming the got set and calculating average of the remaining elements. So, let us have a look at 2D alpha-trimmed mean filter programming. For 2D case we choose window of 3×3 size.

//   2D ALPHA-TRIMMED MEAN FILTER implementation
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
//     alpha  - filter alpha parameter
void _alphatrimmedmeanfilter(const element* image, element* result,
   int N, int M, int alpha)
{
   //   Start of the trimmed ordered set
   const int start = alpha >> 1;
   //   End of the trimmed ordered set
   const int end = 9 - (alpha >> 1);
   //   Move window through all elements of the image
   for (int m = 1; m < M - 1; ++m)
      for (int n = 1; n < N - 1; ++n)
      {
         //   Pick up window elements
         int k = 0;
         element window[9];
         for (int j = m - 1; j < m + 2; ++j)
            for (int i = n - 1; i < n + 2; ++i)
               window[k++] = image[j * N + i];
         //   Order elements (only necessary part of them)
         for (int j = 0; j < end; ++j)
         {
            //   Find position of minimum element
            int min = j;
            for (int l = j + 1; l < 9; ++l)
            if (window[l] < window[min])
               min = l;
            //   Put found minimum element in its place
            const element temp = window[j];
            window[j] = window[min];
            window[min] = temp;
         }
         //   Target index in result image
         const int target = (m - 1) * (N - 2) + n - 1;
         //   Get result - the mean value of the elements in trimmed set
         result[target] = window[start];
         for (int j = start + 1; j < end; ++j)
            result[target] += window[j];
         result[target] /= 9 - alpha;
      }
}

6. Treating edges in 2D case

As for 1D case in 2D case we should extend our input image as well. To do that we are to add lines at the top and at the bottom of the image and add columns to the left and to the right.

Fig. 3. Image extension.Fig. 3. Image extension.

Here is our wrapper function, that does that job.

//   2D ALPHA-TRIMMED MEAN FILTER wrapper
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
//     alpha  - filter alpha parameter
void alphatrimmedmeanfilter(element* image, element* result,
   int N, int M, int alpha)
{
   //   Check arguments
   if (!image || N < 1 || M < 1 || alpha < 0 || 8 < alpha || alpha & 1)
      return;
   //   Allocate memory for signal extension
   element* extension = new element[(N + 2) * (M + 2)];
   //   Check memory allocation
   if (!extension)
      return;
   //   Create image extension
   for (int i = 0; i < M; ++i)
   {
      memcpy(extension + (N + 2) * (i + 1) + 1,
         image + N * i,
         N * sizeof(element));
      extension[(N + 2) * (i + 1)] = image[N * i];
      extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1];
   }
   //   Fill first line of image extension
   memcpy(extension, extension + N + 2, (N + 2) * sizeof(element));
   //   Fill last line of image extension
   memcpy(extension + (N + 2) * (M + 1),
      extension + (N + 2) * M,
      (N + 2) * sizeof(element));
   //   Call alpha-trimmed mean filter implementation
   _alphatrimmedmeanfilter(extension, result ? result : image,
      N + 2, M + 2, alpha);
   //   Free memory
   delete[] extension;
}

7. Alpha-trimmed mean filter library

Now we have four functions, two of them are for processing 1D signals by alpha-trimmed mean filter, and other two are for filtering 2D images. It is time to put everything together and create small alpha-trimmed mean filter library. Let us write code for header file.

//   alphatrimmedmeanfilter.h - declarations for 
//   1D and 2D alpha-trimmed mean filter routines
//
//   The code is property of LIBROW
//   You can use it on your own
//   When utilizing credit LIBROW site

#ifndef _ALPHATRIMMEDMEANFILTER_H_
#define _ALPHATRIMMEDMEANFILTER_H_

//   Signal/image element type
typedef double element;

//   1D ALPHA-TRIMMED MEAN FILTER, window size 5
//     signal - input signal
//     result - output signal, NULL for inplace processing
//     N      - length of the signal
//     alpha  - filter alpha parameter
void alphatrimmedmeanfilter(element* signal, element* result,
   int N, int alpha);

//   2D ALPHA-TRIMMED MEAN FILTER, window size 3x3
//     image  - input image
//     result - output image, NULL for inplace processing
//     N      - width of the image
//     M      - height of the image
//     alpha  - filter alpha parameter
void alphatrimmedmeanfilter(element* image, element* result,
   int N, int M, int alpha);

#endif

Our library is ready. The code we have written is good both for Linux/Unix and Windows platforms. You can download full alpha-trimmed mean filter library source code here:

Download alpha-trimmed mean filter C++ source code (zip, 2 Kb)

Full listings of library files are available online as well:

And now — a couple of applications to play around!

8. Alpha-trimmed mean filter: image restoration

Download alpha-trimmed mean filter for Win32 (zip, 603 Kb)

We have created a couple of applications to show alpha-trimmed mean filter capabilities in restoration images corrupted by impulse noise. The demo package includes 4 files — two applications, sample image and description:

  • alphatrimmedmean.exe — alpha-trimmed mean filter,
  • corrupter.exe — impulse noise generator,
  • sample.bmp — 24-bit sample image,
  • readme.txt — description.

Be aware of the fact, that this sample uses OpenGL, so it should be supported by your system (usually that is the case).

9. Step 1: prepare corrupted image

We have created impulse noise generator that will help us to prepare corrupted image. Start up corrupter.exe application and load image to be corrupted. Choose Set >> Corruption… or click N button in toolbar and set noise level at 5–15%. Click OK. Then save corrupted image.

Fig. 4. Corruption by impulse noise. Screenshot.Fig. 4. Corruption by impulse noise.

10. Step 2: restore corrupted image

Start up alphatrimmedmean.exe application. Load the saved corrupted image. Set alpha parameter to 6: Set >> Alpha… or click α button in toolbar and select 6 in dialog’s combo-box, click OK.

Fig. 5. Alpha parameter selection. Screenshot.Fig. 5. Alpha parameter selection.

Apply alpha-trimmed mean filter by choosing Set >> Filter or clicking F-button in toolbar. See the result. If necessary, filter the image once more.

Fig. 6. Image restored by alpha-trimmed mean filter. Screenshot.Fig. 6. Image restored by alpha-trimmed mean filter.

Posted in Image / Video Filters | Tagged: | Leave a Comment »

Filter window, or filter mask

Posted by Hemprasad Y. Badgujar on October 12, 2014


Filter window, or filter mask

Category. Digital signal and image processing (DSP and DIP) software development.

Abstract. The article deals with the basic concept of signal and image processing science — filter window, or filter mask. Article explains what filter window is and demonstrates window idea for 1D, 2D and 3D cases.

Reference. Have a look at practical utilization of window idea for simple case: Mean filter, or average filter.

3D filter window or mask

1. Introduction to filter window, or filter mask

Practically filter window or filter mask is just a selection pattern. Usually digital filter performs the same actions for every element of the signal or image. For its actions filter uses the neighbor elements as input. Every element has its own neighborhood but the pattern for picking up elements around is the same.

Simple filters usually have one fixed window. More complicated filters, for instance, adaptive ones, can switch window size on fly according to the prior estimation they make. Other filters, for example, hybrid ones, can apply several selection patterns to the same neighborhood — they use several masks.

Now let us look at filter window closer.

2. Filter window or mask

Let us imagine, you should read a letter and what you see in text restricted by hole in a special stencil like this.

First stencilFig. 1. First stencil.

So, the result of reading is sound [t]. Ok, let us read the letter again, but with the help of another stencil:

Second stencilFig. 2. Second stencil.

Now the result of reading t is sound [ð]. Let us make the third try:

Third stencilFig. 3. Third stencil.

Now you are reading letter t as sound [θ].

What happens here? To say that in mathematical language, you are making an operation (reading) over element (letter t). And the result (sound) depends on the element neighborhood (letters next to t).

And that stencil, which helps to pick up element neighborhood, is window! Yes, window is just a stencil or pattern, by means of which you are selecting the element neighborhood — a set of elements around the given one — to help you make decision. Another name for filter window is mask — mask is a stencil, which hides elements we are not paying attention to.

In our example the element we are operating on positioned at leftmost of the window, in practice however its usual position is the center of the window.

Let us see some window examples. In one dimension.

Window or mask of size 5 in 1DFig. 4. Window or mask of size 5 in 1D.

In two dimensions.

Window or mask of size 3x3 in 2DFig. 5. Window or mask of size 3×3 in 2D.

In three dimensions… Think about building. And now — about room in that building. The room is like 3D window, which cuts out some subspace from the entire space of the building. You can find 3D window in volume (voxel) image processing.

Window of mask of size 3x3x3 in 3DFig. 6. Window or mask of size 3×3×3 in 3D.

All our sample windows are of square form, and in practice digital filters usually utilizing that form. But indeed the window form could be any. For example, image filter can use x-like mask for picking up elements.

 

Posted in Image / Video Filters | Tagged: , | Leave a Comment »

 
Extracts from a Personal Diary

dedicated to the life of a silent girl who eventually learnt to open up

Num3ri v 2.0

I miei numeri - seconda versione

ThuyDX

Just another WordPress.com site

Algunos Intereses de Abraham Zamudio Chauca

Matematica, Linux , Programacion Serial , Programacion Paralela (CPU - GPU) , Cluster de Computadores , Software Cientifico

josephdung

thoughts...

Tech_Raj

A great WordPress.com site

Travel tips

Travel tips

Experience the real life.....!!!

Shurwaat achi honi chahiye ...

Ronzii's Blog

Just your average geek's blog

Karan Jitendra Thakkar

Everything I think. Everything I do. Right here.

VentureBeat

News About Tech, Money and Innovation

Chetan Solanki

Helpful to u, if u need it.....

ScreenCrush

Explorer of Research #HEMBAD

managedCUDA

Explorer of Research #HEMBAD

siddheshsathe

A great WordPress.com site

Ari's

This is My Space so Dont Mess With IT !!

%d bloggers like this: