// This program is free software: you can use, modify and/or redistribute it
// under the terms of the simplified BSD License. You should have received a
// copy of this license along this program. If not, see
// <http://www.opensource.org/licenses/bsd-license.html>.
//
// Copyright (C) 2012, Javier Sánchez Pérez <jsanchez@dis.ulpgc.es>
// Copyright (C) 2014, Nelson Monzón López <nmonzon@ctim.es>
// All rights reserved.
#ifndef MASK_H
#define MASK_H
#include <omp.h>
/**
*
* Function to apply a 3x3 mask to an image
*
*/
void
mask3x3 (const float *input, //input image
float *output, //output image
const int nx, //image width
const int ny, //image height
const int nz, // number of color channels in the image
const float *mask //mask to be applied
)
{
int nx_multichannel = nx * nz;
for(int index_multichannel = 0; index_multichannel < nz; index_multichannel++){
//apply the mask to the center body of the image
#pragma omp parallel for
for (int i = 1; i < ny - 1; i++){
for (int j = 1; j < nx - 1; j++){
int k = (i * nx + j) * nz + index_multichannel;
double sum = 0;
for (int l = 0; l < 3; l++){
for (int m = 0; m < 3; m++){
int p = ((i + l - 1) * nx + j + m - 1) * nz + index_multichannel;
sum += input[p] * mask[l * 3 + m];
}
}
output[k] = sum;
}
}
//apply the mask to the first and last rows
#pragma omp parallel for
for (int j = 1; j < nx - 1; j++) {
int index = j * nz + index_multichannel;
double sum = 0;
sum += input[index - nz] * (mask[0] + mask[3]);
sum += input[index] * (mask[1] + mask[4]);
sum += input[index + nz] * (mask[2] + mask[5]);
sum += input[nx_multichannel + j - nz] * mask[6];
sum += input[nx_multichannel + j] * mask[7];
sum += input[nx_multichannel + j + nz] * mask[8];
output[j] = sum;
index = ((ny - 2) * nx + j) * nz + index_multichannel;
sum = 0;
sum += input[index - nz] * mask[0];
sum += input[index] * mask[1];
sum += input[index + nz] * mask[2];
index = ((ny - 1) * nx + j) * nz + index_multichannel;
sum += input[index - nz] * (mask[6] + mask[3]);
sum += input[index] * (mask[7] + mask[4]);
sum += input[index + 1] * (mask[8] + mask[5]);
output[index] = sum;
}
//apply the mask to the first and last columns
#pragma omp parallel for
for (int i = 1; i < ny - 1; i++){
int index = i * nx_multichannel + index_multichannel;
double sum = 0;
int index_row = (i - 1) * nx_multichannel + index_multichannel;
sum += input[index_row] * (mask[0] + mask[1]);
sum += input[index_row + nz] * mask[2];
sum += input[index] * (mask[3] + mask[4]);
sum += input[index + nz] * mask[5];
index_row = (i + 1) * nx_multichannel + index_multichannel;
sum += input[index_row] * (mask[6] + mask[7]);
sum += input[index_row + nz] * mask[8];
output[index] = sum;
sum = 0;
sum += input[index - 2 * nz] * mask[0];
sum += input[index - nz] * (mask[1] + mask[2]);
index_row = (i + 1) * nx_multichannel + index_multichannel;
sum += input[index_row - 2 * nz] * mask[3];
sum += input[index_row - nz] * (mask[4] + mask[5]);
index_row = (i + 2) * nx_multichannel + index_multichannel;
sum += input[index_row - 2 * nz] * mask[6];
sum += input[index_row - nz] * (mask[7] + mask[8]);
output[(i * nx + nx - 1) * nz + index_multichannel] = sum;
}
//apply the mask to the four corners
output[index_multichannel] =
input[index_multichannel] * (mask[0] + mask[1] + mask[3] + mask[4]) +
input[index_multichannel + nz] * (mask[2] + mask[5]) +
input[nx_multichannel + index_multichannel] * (mask[6] + mask[7]) +
input[nx_multichannel + index_multichannel + nz] * mask[8];
output[nx_multichannel - nz + index_multichannel] =
input[(nx - 2) * nz + index_multichannel] * (mask[0] + mask[3]) +
input[(nx - 1) * nz + index_multichannel] * (mask[1] + mask[2] + mask[4] + mask[5]) +
input[(2 * nx - 2) * nz + index_multichannel] * mask[6] +
input[(2 * nx - 1) * nz + index_multichannel] * (mask[7] + mask[8]);
output[(ny - 1) * nx_multichannel + index_multichannel] =
input[(ny - 2) * nx_multichannel + index_multichannel] * (mask[0] + mask[1]) +
input[((ny - 2) * nx + 1) * nz + index_multichannel] * mask[2] +
input[(ny - 1) * nx_multichannel + index_multichannel] * (mask[3] + mask[4] + mask[6] + mask[7]) +
input[((ny - 1) * nx + 1) * nz + index_multichannel] * (mask[5] + mask[8]);
output[(ny * nx - 1) * nz + index_multichannel] =
input[((ny - 1) * nx - 2) * nz + index_multichannel] * mask[0] +
input[((ny - 1) * nx - 1) * nz + index_multichannel] * (mask[1] + mask[2]) +
input[(ny * nx - 2) * nz + index_multichannel] * (mask[3] + mask[6]) +
input[(ny * nx - 1) * nz + index_multichannel] * (mask[4] + mask[5] + mask[7] + mask[8]);
} // end loop for channels information
} // end mask3x3
/**
*
* Compute the second order X derivative
*
*/
void
Dxx (const float *I, //input image
float *Ixx, //oputput derivative
const int nx, //image width
const int ny, //image height
const int nz //number of color channels in the image
)
{
//mask of second derivative
float M[] = { 0., 0., 0.,
1., -2., 1.,
0., 0., 0.
};
//computing the second derivative
mask3x3 (I, Ixx, nx, ny, nz, M);
}
/**
*
* Compute the second order Y derivative
*
*/
void
Dyy (const float *I, //input image
float *Iyy, //oputput derivative
const int nx, //image width
const int ny, //image height
const int nz //number of color channels in the image
)
{
//mask of second derivative
float M[] = { 0., 1., 0.,
0., -2., 0.,
0., 1., 0.
};
//computing the second derivative
mask3x3 (I, Iyy, nx, ny, nz, M);
}
/**
*
* Compute the second order XY derivative
*
*/
void
Dxy (const float *I, //input image
float *Ixy, //oputput derivative
const int nx, //image width
const int ny, //image height
const int nz //number of color channels in the image
)
{
//mask of second derivative
float M[] = { 1. / 4., 0., -1. / 4.,
0., 0., 0.,
-1. / 4., 0., 1. / 4.
};
//computing the second derivative
mask3x3 (I, Ixy, nx, ny, nz, M);
}
/**
*
* Compute the gradient with central differences
*
*/
void
gradient (const float *input, //input image
float *dx, //computed x derivative
float *dy, //computed y derivative
const int nx, //image width
const int ny, //image height
const int nz //number of color channels in the image
)
{
const int nx_multichannel = nx * nz;
for(int index_multichannel = 0; index_multichannel < nz; index_multichannel++){
//gradient in the center body of the image
#pragma omp parallel for
for (int i = 1; i < ny - 1; i++){
for (int j = 1; j < nx - 1; j++){
const int k = (i * nx + j) * nz + index_multichannel;
dx[k] = 0.5 * (input[k + nz] - input[k - nz]);
dy[k] = 0.5 * (input[k + nx_multichannel] - input[k - nx_multichannel]);
}
}
//gradient in the first and last rows
#pragma omp parallel for
for (int j = 1; j < nx - 1; j++){
const int index = j * nz + index_multichannel;
dx[j] = 0.5 * (input[index + nz] - input[index - nz]);
dy[j] = 0.5 * (input[index + nx_multichannel] - input[index]);
const int k = ((ny - 1) * nx + j) * nz + index_multichannel;
dx[k] = 0.5 * (input[k + nz] - input[k - nz]);
dy[k] = 0.5 * (input[k] - input[k - nx_multichannel]);
}
//gradient in the first and last columns
#pragma omp parallel for
for (int i = 1; i < ny - 1; i++) {
const int p = (i * nx_multichannel) + index_multichannel;
dx[p] = 0.5 * (input[p + nz] - input[p]);
dy[p] = 0.5 * (input[p + nx_multichannel] - input[p - nx_multichannel]);
const int k = ((i + 1) * nx - 1) * nz + index_multichannel;
dx[k] = 0.5 * (input[k] - input[k - nz]);
dy[k] = 0.5 * (input[k + nx_multichannel] - input[k - nx_multichannel]);
}
//calculate the gradient in the corners
dx[index_multichannel] = 0.5 * (input[index_multichannel + nz] - input[index_multichannel]);
dy[index_multichannel] = 0.5 * (input[nx_multichannel + index_multichannel] - input[index_multichannel]);
const int corner_up_right = (nx-1) * nz + index_multichannel;
dx[corner_up_right] = 0.5 * (input[corner_up_right] - input[corner_up_right - nz]);
dy[corner_up_right] = 0.5 * (input[(2 * nx_multichannel) + index_multichannel - nz] - input[corner_up_right]);
const int corner_down_left = ((ny - 1) * nx) * nz + index_multichannel;
dx[corner_down_left] = 0.5 * (input[corner_down_left + nz] - input[corner_down_left]);
dy[corner_down_left] = 0.5 * (input[corner_down_left] - input[(ny - 2) * nx_multichannel + index_multichannel]) ;
const int corner_down_right = ny * nx_multichannel - nz + index_multichannel;
dx[corner_down_right] = 0.5 * (input[corner_down_right] - input[corner_down_right - nz]);
dy[corner_down_right] = 0.5 * (input[corner_down_right] - input[(ny - 1) * nx_multichannel - nz + index_multichannel]);
} // end loop for multi-channel
} // end gradient function
#endif