Deconvolution

// freeimage_test2014.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <FreeImage.h>
#include <string>
#pragma comment(lib,"FreeImage.lib")
#pragma comment(lib,"FreeImage.dll")

#define IN
#define OUT
#define UNKNOWN 0
#define BMP 1
#define GIF 2
#define JPEG 3
#define PNG 4
#define TIFF 5
#define BYTE unsigned char
#define FOUR_SIDES 1
#define EIGHT_SIDES 2

struct IMAGE {
	BYTE *red; // Пиксельная карта красного канала
	BYTE *green; // Пиксельная карта зеленого канала
	BYTE *blue; // Пиксельная карта синего канала
	BYTE *lum; // Пиксельная карта полутонового изображения (в градациях серого)
	int width; // Ширина изображения
	int height; // Высота изоражения
};

IMAGE loadImage(const char *name, int type) {
	IMAGE result;
	FIBITMAP *bitmap;
	FREE_IMAGE_FORMAT fif;

	result.blue = 0;
	result.height = 0;
	result.width = 0; 
	result.lum = 0;
	result.red = 0;
	result.green = 0;

	switch(type) {
		case 1: 
			fif = FIF_BMP; 
			break;
		case 2: 
			fif = FIF_GIF; 
			break;
		case 3: 
			fif = FIF_JPEG; 
			break;
		case 4: 
			fif = FIF_PNG; 
			break;
		case 5: 
			fif = FIF_TIFF; 
			break;
		default:
			fif = FIF_UNKNOWN;
	}
	if (fif == FIF_UNKNOWN) {
		printf("loadImage: unknown file format: %d\n", type);
	} else {
		bitmap = FreeImage_Load(fif, name, 0); // Загружаем изображение в память
		if (bitmap) {
			printf("loadImage: image was loaded! Year!\n");
			int width, height; // Ширина и высота загруженного изображения в пикселах
			BYTE *red_map; // Указатель на красную пиксельную карту
			BYTE *green_map; // Указатель на зеленую пиксельную карту
			BYTE *blue_map; // Указатель на синюю пиксельную карту
			int x, y; // Координаты пикселя
			bool success; // Удалось ли получить все пиксели
			RGBQUAD color; // Цвет пикселя
			
			width = FreeImage_GetWidth(bitmap);
			height = FreeImage_GetHeight(bitmap);
			
			result.red = new BYTE[width*height];
			result.green = new BYTE[width*height];
			result.blue = new BYTE[width*height];
			result.width = width;
			result.height = height;
			
			success = true;
			for (x = 0; x < width; x++) {
				for (y = 0; y < height; y++) {
					success = FreeImage_GetPixelColor(bitmap, x, y, &color);
					result.red[width*y + x] = color.rgbRed;
					result.green[width*y + x] = color.rgbGreen;
					result.blue[width*y + x] = color.rgbBlue;
					if (!success) break;
				}
				if (!success) break;
			}
			if (success) {
				printf("loadImage: color maps were formed!\n");
			} else {
				printf("loadImage: couldn\'t get pixel (%d, %d)\n", x, y);
			}
			FreeImage_Unload(bitmap);
		} else {
			printf("loadImage: image was not loaded\n");
		}
	}
	return result;
}

void saveImage(IMAGE image, const char *name, int type) {
	FIBITMAP *bitmap;
	FREE_IMAGE_FORMAT fif;
	switch(type) {
		case 1: 
			fif = FIF_BMP; 
			break;
		case 2: 
			fif = FIF_GIF; 
			break;
		case 3: 
			fif = FIF_JPEG; 
			break;
		case 4: 
			fif = FIF_PNG; 
			break;
		case 5: 
			fif = FIF_TIFF; 
			break;
		default:
			fif = FIF_UNKNOWN;
	}
	if (fif == FIF_UNKNOWN) {
		printf("saveImage: unknown file format: %d\n", type);
		return;
	}
	bitmap = FreeImage_Allocate(image.width, image.height, 24); // Выделяем память под изображение
	if (bitmap) {
		printf("saveImage: image was created!\n");
		int width, height; // Ширина и высота загруженного изображения в пикселах
		int x, y; // Координаты пикселя
		bool success; // Удалось ли получить все пиксели
		RGBQUAD color; // Цвет пикселя
		
		width = image.width;
		height = image.height;

		success = true;
		for (x = 0; x < width; x++) {
			for (y = 0; y < height; y++) {
				color.rgbRed = image.red[width*y + x];
				color.rgbGreen = image.green[width*y + x];
				color.rgbBlue = image.blue[width*y + x];
				success = FreeImage_SetPixelColor(bitmap, x, y, &color);
				if (!success) break;
			}
			if (!success) break;
		}
		if (success) {
			printf("saveImage: color maps were writed into bitmap!\n");
			if (FreeImage_Save(fif, bitmap, name)) {
				printf("saveImage: bitmap is successfully saved!\n");
			} else {
				printf("saveImage: bitmap couldn\'t be saved\n");
			}
		} else {
			printf("saveImage: couldn\'t set pixel (%d, %d)\n", x, y);
		}
		FreeImage_Unload(bitmap);
	} else {
		printf("saveImage: image was not created\n");
	}
}
//-------------------------------------------------

void toGrayscale(IMAGE *image) {
	int w = image->width;
	int h = image->height;
	int i, lum;
	for (i = 0; i < w*h; i++) {
		lum = 0.299*image->red[i] + 0.587*image->green[i] + 0.114*image->blue[i];
		if (lum > 255) lum = 255;
		image->red[i] = lum;
	}
	delete [] image->green;
	delete [] image->blue;
	image->lum = image->red;
	image->green = image->red;
	image->blue = image->red;
}
//-------------------------------------------------


void inverse(IMAGE image) {
	int i;
	for (i = 0; i < image.width*image.height; i++) {
		image.red[i] = 255 - image.red[i];
		image.green[i] = 255 - image.green[i];
		image.blue[i] = 255 - image.blue[i];
	}
}
//-------------------------------------------------

IMAGE *conv(IMAGE image, IMAGE psf) {
	if (!psf.lum) {
		printf("conv: PSF should be a grayscale image\n");
		return 0;
	}
	if (psf.width%2 != 0 || psf.height%2 != 0) {
		printf("conv: PSF should have a center\n");
		return 0;
	}
	int w1 = image.width, h1 = image.height;
	int w2 = psf.width, h2 = psf.height;
	int a = psf.width/2; // Горизонтальное расстояние от центра psf до края
	int b = psf.height/2; // Вертикальное расстояние от центра psf до края
	int x, y, i, j;
	for (x = 0; x < w1; x++) {
		for (y = 0; y < h1; y++) { // Проход по каждой точке изображения f
			double sum = 0;
			for (i = 0; i < w2; i++) {
				for (j = 0; j < h2; j++) { // Проход по каждной точке ФРТ g
					int index_x = x + i - a;
					int index_y = y + j - b;
					if (index_x < 0) index_x = 0; if (index_x >= w1) index_x = w1 - 1;
					if (index_y < 0) index_y = 0; if (index_y >= h1) index_y = h1 - 1;
					sum += g[i + j*w2]*f[index_x + index_y*w1];
				}
			}
			dest[x + y*w1] = sum;
		}
	}
}

void convolution(double *f, int w1, int h1, double *g, int w2, int h2, double *dest) {
	int a = ceil((double)(w2 - 1)/2); // Горизонтальное расстояние от центра psf до края
	int b = ceil((double)(h2 - 1)/2); // Вертикальное расстояние от центра psf до края
	for (int x = 0; x < w1; x++) {
		for (int y = 0; y < h1; y++) { // Проход по каждой точке изображения f
			double sum = 0;
			
			for (int i = 0; i < w2; i++) {
				for (int j = 0; j < h2; j++) { // Проход по каждной точки маски g
					int index_x = x + i - a;
					int index_y = y + j - b;
					if (index_x < 0) index_x = 0; if (index_x >= w1) index_x = w1 - 1;
					if (index_y < 0) index_y = 0; if (index_y >= h1) index_y = h1 - 1;
					sum += g[i + j*w2]*f[index_x + index_y*w1];
				}
			}
			dest[x + y*w1] = sum;
		}
	}
}
//--------------------------

void rich(double *g, int w1, int h1, double *h, int w2, int h2, double *dest, int iterations) {
	double *prev = new double[w1*h1];
	double *temp = new double[w1*h1];
	for (int i = 0; i < w1*h1; i++) prev[i] = dest[i] = g[i];
	for (int k = 0; k < iterations; k++) {
		convolution(prev, w1, h1, h, w2, h2, temp);
		for (int i = 0; i < w1*h1; i++) if (temp[i]) temp[i] = g[i]/temp[i];
		convolution(temp, w1, h1, h, w2, h2, dest);
		for (int i = 0; i < w1*h1; i++) prev[i] = prev[i]*dest[i];
	}
	for (int i = 0; i < w1*h1; i++) dest[i] = prev[i];

	delete [] prev;
	delete [] temp;
}
//---------------------------------------------------

void laplace(IMAGE image, int type) {
	BYTE *map, *buf;
	int j, i, x, y, w = image.width, h = image.height;

	for (j = 0; j < 3; j++) { // Цикл по цветовым каналам
		switch(j) {
			case 0: map = image.red; printf("red\n"); break;
			case 1: map = image.green; printf("green\n"); break;
			case 2: map = image.blue; printf("blue\n"); break;
		}

		buf = new BYTE[w*h];
		for (i = 0; i < w*h; i++) {
			buf[i] = map[i];
		}
		
		for (x = 1; x < w-1; x++) {
			for (y = 1; y < h-1; y++) {
				if (type == (type|FOUR_SIDES)) {
					int lum = 5*buf[y*w+x];
					lum -= buf[y*w+x+1] + buf[y*w+x-1] + buf[(y+1)*w+x] + buf[(y-1)*w+x];
					if (lum < 0) lum = 0;
					if (lum > 255) lum = 255;
					map[y*w+x] = lum;
				} 
				else {
					int lum = 9*buf[y*w+x];
					lum -= buf[y*w+x+1] + buf[y*w+x-1] + buf[(y+1)*w+x] + buf[(y-1)*w+x];
					lum -= buf[(y+1)*w+x+1] + buf[(y+1)*w+x-1] + buf[(y-1)*w+x+1] + buf[(y-1)*w+x-1];
					if (lum < 0) lum = 0;
					if (lum > 255) lum = 255;
					map[y*w+x] = lum;
				}
			}
		}
		if (image.lum) {
			printf("laplace: image is grayscale\n");
			break; // Если изображение полутоновое, не переходить к другим цветовым каналам
		}
	}
}
//--------------------------------------------

int main()
{
	IMAGE image;
	BYTE *map;
	image = loadImage("images/no_noise.png", PNG);
	
	//toGrayscale(&image);
	laplace(image, FOUR_SIDES);

	saveImage(image, "images/no_noise_res.png", PNG);
	return 0;
}