1 | initial version |
I have not used the OpenCV to implement it.I've already implemented this code in ANSI C language. This code simply by changing the name of its functions is to become the Opencv C Edition.Stationary Gaussian kernel is used in the regular version but in the gaussian filter with mask for each pixel non zero in the mask you have to compute gaussian kernel for it.
typedef struct BhImageA
{
int width; //Image width in pixels
int height; // Image height in pixels
int widthStep; // The size of an aligned image row, in bytes
byte* imageData; //A pointer to the aligned image data
byte depth; //Channel depth in bits
BhRectA roi; // Region Of Interest (ROI). If not NULL, only this image region will be processed
}BhImageA;
typedef struct BhRectA
{
int x; // Offset (usually the top-left corner) and size of a rectangle
int y; // y-coordinate of the top-left corner
int width; // y-coordinate of the top-left corner
int height; // y-coordinate of the top-left corner
}BhRectA;
#define BH_DEPTH_8 8
#define BH_DEPTH_32 32
#define BH_DEPTH_32_F 64
BhSizeA bhSizeA(int width,int height)
{
BhSizeA result;
result.width = width;
result.height = height;
return result;
}
BhRectA bhGetImageRectA(BhImageA* srcImage)
{
BhRectA result;
result.x = 0;
result.y = 0;
result.width = srcImage->width;
result.height = srcImage->height;
return result;
}
int bhGetWidthStepA(int width)
{
return int(ceil(width / 4.0)*4);
}
void bhMulSA(BhImageA* srcImage,BhImageA* dstImage,double value)
{
if (srcImage->depth == BH_DEPTH_8)
{
byte v=(byte)value;
for (int i = 0; i < srcImage->height; i++)
{
byte* srcRow = srcImage->imageData + srcImage->widthStep * i;
byte* dstRow = dstImage->imageData + dstImage->widthStep * i;
for (int j = 0; j < srcImage->width; j++)
dstRow[j] = srcRow[j]* v;
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
int v=(int)value;
for (int i = 0; i < srcImage->height; i++)
{
int* srcRow = (int*)srcImage->imageData + srcImage->widthStep * i;
int* dstRow = (int*)dstImage->imageData + dstImage->widthStep * i;
for (int j = 0; j < srcImage->width; j++)
dstRow[j] = srcRow[j]* v;
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float v=(double)value;
for (int i = 0; i < srcImage->height; i++)
{
float* srcRow = (float*)srcImage->imageData + srcImage->widthStep * i;
float* dstRow = (float*)dstImage->imageData + dstImage->widthStep * i;
for (int j = 0; j < srcImage->width; j++)
dstRow[j] = srcRow[j]* v;
}
}
}
double bhSumA(BhImageA* srcImage,BhImageA* maskImage)
{
double sum = 0;
BhRectA roi = srcImage->roi;
if (maskImage == NULL)
{
if (srcImage->depth == BH_DEPTH_8)
{
for (int i = roi.y; i < roi.y + roi.height; i++)
{
byte* row = srcImage->imageData + srcImage->widthStep * i;
for (int j = roi.x; j < roi.x + roi.width; j++)
sum += row[j];
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
for (int i = roi.y; i < roi.y + roi.height ; i++)
{
int* row = (int*)srcImage->imageData + srcImage->widthStep * i;
for (int j = roi.x ; j < roi.x + roi.width; j++)
sum += row[j];
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
for (int i = roi.y; i < roi.y + roi.height ; i++)
{
float* row = (float*)srcImage->imageData + srcImage->widthStep * i;
for (int j = roi.x; j < roi.x + roi.width; j++)
sum += row[j];
}
}
}
else if (maskImage != NULL)
{
if (srcImage->width != maskImage->width && srcImage->height != maskImage->height)
return -1;
if (srcImage->depth == BH_DEPTH_8)
{
for (int i = roi.y; i < roi.y + roi.height; i++)
{
byte* row = srcImage->imageData + srcImage->widthStep * i;
byte* rowMask = maskImage->imageData + maskImage->widthStep * i;
for (int j = roi.x; j < roi.x + roi.width; j++)
if (rowMask[j] != 0)
sum += row[j];
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
for (int i = roi.y; i < roi.y + roi.height ; i++)
{
int* row = (int*)srcImage->imageData + srcImage->widthStep * i;
byte* rowMask = maskImage->imageData + maskImage->widthStep * i;
for (int j = roi.x ; j < roi.x + roi.width; j++)
if (rowMask[j] != 0)
sum += row[j];
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
for (int i = roi.y; i < roi.y + roi.height ; i++)
{
float* row = (float*)srcImage->imageData + srcImage->widthStep * i;
byte* rowMask = maskImage->imageData + maskImage->widthStep * i;
for (int j = roi.x; j < roi.x + roi.width; j++)
if (rowMask[j] != 0)
sum += row[j];
}
}
}
return sum;
}
void bhSetImageROIA(BhImageA* srcImage,BhRectA roi)
{
srcImage->roi = roi;
}
void bhResetImageROIA(BhImageA* srcImage)
{
srcImage->roi = bhGetImageRectA(srcImage);
}
int bhCopyImageA(BhImageA* srcImage , BhImageA* dstImage)
{
if ( (srcImage->depth != dstImage->depth) || (srcImage->roi.width != dstImage->roi.width) || (srcImage->roi.height != dstImage->roi.height))
return -1;
if (srcImage->depth == BH_DEPTH_8)
{
byte* srcImageData = srcImage->imageData;
byte* srcRow = NULL;
byte* dstImageData = dstImage->imageData;
byte* dstRow = NULL;
BhRectA srcRoi = srcImage->roi;
BhRectA dstRoi = dstImage->roi;
for (int i = 0 ; i < srcRoi.height ;i++)
{
srcRow = srcImageData + srcImage->widthStep * (i + srcRoi.y) ;
dstRow = dstImageData + dstImage->widthStep * (i + dstRoi.y) ;
for (int j=0; j < srcRoi.width; j++)
{
dstRow[j + dstRoi.x] = srcRow[j + srcRoi.x];
}
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
int* srcImageData = (int*)srcImage->imageData;
int* srcRow = NULL;
int* dstImageData = (int*)dstImage->imageData;
int* dstRow = NULL;
BhRectA srcRoi = srcImage->roi;
BhRectA dstRoi = dstImage->roi;
for (int i = 0 ; i < srcRoi.height ;i++)
{
srcRow = srcImageData + srcImage->widthStep * (i + srcRoi.y) ;
dstRow = dstImageData + dstImage->widthStep * (i + dstRoi.y) ;
for (int j=0; j < srcRoi.width; j++)
{
dstRow[j + dstRoi.x] = srcRow[j + srcRoi.x];
}
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float* srcImageData = (float*)srcImage->imageData;
float* srcRow = NULL;
float* dstImageData = (float*)dstImage->imageData;
float* dstRow = NULL;
BhRectA srcRoi = srcImage->roi;
BhRectA dstRoi = dstImage->roi;
for (int i = 0 ; i < srcRoi.height ;i++)
{
srcRow = srcImageData + srcImage->widthStep * (i + srcRoi.y) ;
dstRow = dstImageData + dstImage->widthStep * (i + dstRoi.y) ;
for (int j=0; j < srcRoi.width; j++)
{
dstRow[j + dstRoi.x] = srcRow[j + srcRoi.x];
}
}
}
return 0;
}
BhImageA* bhCreateImageA(BhSizeA size,byte depth)
{
BhImageA* result = (BhImageA*)malloc(sizeof(BhImageA));
result->width = size.width;
result->height = size.height;
result->depth = depth;
if (depth == BH_DEPTH_8)
{
result->widthStep = bhGetWidthStepA(size.width);
result->imageData = (byte*)malloc(result->widthStep * size.height);
}
else if (depth == BH_DEPTH_32)
{
result->widthStep = size.width;
result->imageData = (byte*)malloc(size.width * sizeof(int) * size.height);
}
else if (depth == BH_DEPTH_32_F)
{
result->widthStep = size.width;
result->imageData = (byte*)malloc(size.width * sizeof(float) * size.height);
}
result->roi = bhGetImageRectA(result);
return result;
}
void bhReleaseImageA(BhImageA* srcImage)
{
free( srcImage->imageData);
free(srcImage);
}
BhImageA* bhMulTransposedA(BhImageA* srcImage)
{
BhImageA* result = NULL;
if (srcImage->depth == BH_DEPTH_32_F)
{
result = bhCreateImageA(bhSizeA(srcImage->width,srcImage->width),BH_DEPTH_32_F);
float* srcRow = (float*)srcImage->imageData;
for (int i = 0; i < result->height ; i++)
{
float* dstRow = (float*)result->imageData + result->width * i;
for (int j = 0; j < result->width; j++)
{
dstRow[j] = srcRow[i] * srcRow[j] * result->width;
}
}
}
return result;
}
BhImageA* bhGetGaussianKernelA( int n, double sigma )
{
const int SMALL_GAUSSIAN_SIZE = 7;
static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
{
{1.f},
{0.25f, 0.5f, 0.25f},
{0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
{0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
};
const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
small_gaussian_tab[n>>1] : 0;
//CV_Assert( ktype == CV_32F || ktype == CV_64F );
BhImageA* kernel = bhCreateImageA(bhSizeA(n,1),BH_DEPTH_32_F);
float* cf = (float*)kernel->imageData;
double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
double scale2X = -0.5/(sigmaX*sigmaX);
double sum = 0;
int i;
for( i = 0; i < n; i++ )
{
double x = i - (n-1)*0.5;
double t = fixed_kernel ? (double)fixed_kernel[i] : exp(scale2X*x*x);
cf[i] = (float)t;
sum += cf[i];
}
sum = 1./sum;
for( i = 0; i < n; i++ )
{
cf[i] = (float)(cf[i]*sum);
}
return kernel;
}
void bhFillBorderA(BhImageA* srcImage,BhImageA* dstImage,BhSizeA offset,int borderType,float value)
{
if (borderType == BH_BORDER_CONSTANT)
{
if (srcImage->depth == BH_DEPTH_8)
{
byte v = (byte)value;
for (int i = 0; i < offset.height; i++)
{
byte* dstRow1 = dstImage->imageData + dstImage->widthStep * i;
byte* dstRow2 = dstImage->imageData + dstImage->widthStep * (dstImage->height - i - 1);
for (int j = 0; j < dstImage->width; j++)
{
dstRow1[j] = v;
dstRow2[j] = v;
}
}
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
byte* dstRow = dstImage->imageData + dstImage->widthStep * i;
for (int j = 0; j < offset.width; j++)
{
dstRow[j] = v;
dstRow[dstImage->widthStep - j -1] = v;
}
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
int v = (int)value;
for (int i = 0; i < offset.height; i++)
{
int* dstRow1 = (int*)dstImage->imageData + dstImage->width * i;
int* dstRow2 = (int*)dstImage->imageData + dstImage->width * (dstImage->height - i - 1);
for (int j = 0; j < dstImage->width; j++)
{
dstRow1[j] = v;
dstRow2[j] = v;
}
}
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
int* dstRow = (int*)dstImage->imageData + dstImage->width * i;
for (int j = 0; j < offset.width; j++)
{
dstRow[j] = v;
dstRow[dstImage->width - j -1] = v;
}
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float v = value;
for (int i = 0; i < offset.height; i++)
{
float* dstRow1 = (float*)dstImage->imageData + dstImage->width * i;
float* dstRow2 = (float*)dstImage->imageData + dstImage->width * (dstImage->height - i - 1);
for (int j = 0; j < dstImage->width; j++)
{
dstRow1[j] = v;
dstRow2[j] = v;
}
}
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
float* dstRow = (float*)dstImage->imageData + dstImage->width * i;
for (int j = 0; j < offset.width; j++)
{
dstRow[j] = v;
dstRow[dstImage->widthStep - j -1] = v;
}
}
}
}
else if (borderType == BH_BORDER_REPLICATE)
{
if (srcImage->depth == BH_DEPTH_8)
{
byte* srcRow1 = dstImage->imageData + dstImage->widthStep * offset.height + offset.width;
byte* srcRow2 = dstImage->imageData + dstImage->widthStep * (dstImage->height - offset.height-1) + offset.width;
size_t size = srcImage->width;
for (int i = 0; i < offset.height; i++)
{
byte* dstRow1 = dstImage->imageData + dstImage->widthStep * i + offset.width;
memcpy(dstRow1,srcRow1,size);
byte* dstRow2 = dstImage->imageData + dstImage->widthStep * (dstImage->height - i - 1) + offset.width;
memcpy(dstRow2,srcRow2,size);
}
size = offset.width;
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
byte* dstRow = dstImage->imageData + dstImage->widthStep * i;
byte v1 = dstRow[offset.width];
byte v2 = dstRow[dstImage->width - offset.width-1] ;
memset(&dstRow[0],v1,size);
memset(&dstRow[dstImage->width - offset.width],v2,size);
}
}
else if (srcImage->depth == BH_DEPTH_32)
{
int* srcRow1 = (int*)dstImage->imageData + dstImage->width * offset.height + offset.width;
int* srcRow2 = (int*)dstImage->imageData + dstImage->width * (dstImage->height - offset.height-1) + offset.width;
size_t size = srcImage->width * sizeof(int);
for (int i = 0; i < offset.height; i++)
{
int* dstRow1 = (int*)dstImage->imageData + dstImage->width * i + offset.width;
memcpy(dstRow1,srcRow1,size);
int* dstRow2 = (int*)dstImage->imageData + dstImage->width * (dstImage->height - i - 1) + offset.width;
memcpy(dstRow2,srcRow2,size);
}
size = offset.width * sizeof(int);
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
int* dstRow = (int*)dstImage->imageData + dstImage->width * i;
int v1 = dstRow[offset.width];
int v2 = dstRow[dstImage->width - offset.width-1] ;
for (int j = 0; j < offset.width; j++)
{
dstRow[j] = v1;
dstRow[dstImage->width - j] = v1;
}
}
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float* srcRow1 = (float*)dstImage->imageData + dstImage->width * offset.height + offset.width;
float* srcRow2 = (float*)dstImage->imageData + dstImage->width * (dstImage->height - offset.height-1) + offset.width;
size_t size = srcImage->width * sizeof(float);
for (int i = 0; i < offset.height; i++)
{
float* dstRow1 = (float*)dstImage->imageData + dstImage->width * i + offset.width;
memcpy(dstRow1,srcRow1,size);
float* dstRow2 = (float*)dstImage->imageData + dstImage->width * (dstImage->height - i - 1) + offset.width;
memcpy(dstRow2,srcRow2,size);
}
size = offset.width * sizeof(float);
for (int i = offset.height; i < dstImage->height - offset.height; i++)
{
float* dstRow = (float*)dstImage->imageData + dstImage->width * i;
float v1 = dstRow[offset.width];
float v2 = dstRow[dstImage->width - offset.width-1] ;
for (int j = 0; j < offset.width; j++)
{
dstRow[j] = v1;
dstRow[dstImage->width - j] = v1;
}
}
}
}
else if (borderType == BH_BORDER_COPY)
{
int size;
if (srcImage->depth == BH_DEPTH_8)
size = 1;
else if (srcImage->depth == BH_DEPTH_32)
size = sizeof(int);
else if (srcImage->depth = BH_DEPTH_32_F)
size = sizeof(float);
for (int i = 0; i < offset.height; i++)
{
byte* srcRow = srcImage->imageData + srcImage->widthStep * size * i;
byte* dstRow = dstImage->imageData + dstImage->widthStep * size * i;
memcpy(dstRow,srcRow,srcImage->widthStep * size);
srcRow = srcImage->imageData + srcImage->widthStep * size * (srcImage->height - i -1);
dstRow = dstImage->imageData + dstImage->widthStep * size * (srcImage->height - i -1);
memcpy(dstRow,srcRow,srcImage->widthStep * size);
}
for (int i = offset.height; i < srcImage->height - offset.height; i++)
{
byte* srcRow = srcImage->imageData + srcImage->widthStep * size * i;
byte* dstRow = dstImage->imageData + dstImage->widthStep * size * i;
memcpy(dstRow,srcRow,offset.width * size);
memcpy(dstRow + (srcImage->width - offset.width) * size,srcRow+ (srcImage->width - offset.width) * size,offset.width * size);
}
}
}
BhImageA* bhCreateFillBorderA(BhImageA* srcImage,BhSizeA offset,int borderType,float value)
{
BhImageA* resultImg = bhCreateBorderImageA(srcImage,offset);
bhFillBorderA(srcImage,resultImg,offset,borderType,value);
return resultImg;
}
BhImageA* bhCreateBorderImageA(BhImageA* srcImage,BhSizeA borderSize)
{
BhImageA* resultImg = bhCreateImageA(bhSizeA(srcImage->width + borderSize.width * 2,srcImage->height + borderSize.height * 2),srcImage->depth);
bhSetImageROIA(resultImg,bhRectA(borderSize.width,borderSize.height,srcImage->width,srcImage->height));
bhCopyImageA(srcImage,resultImg);
bhResetImageROIA(resultImg);
return resultImg;
}
void** bhCreateRowsA(BhImageA* srcImage)
{
if (srcImage->depth == BH_DEPTH_8)
{
byte** rows = (byte**) malloc(srcImage->height * sizeof(byte*) );
for (int y=0 ;y < srcImage->height;y++)
rows[y] = srcImage->imageData + srcImage->widthStep * y;
return (void**)rows;
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float** rows = (float**) malloc(srcImage->height * sizeof(float*));
for (int y = 0; y < srcImage->height ; y++)
rows[y] = (float*)srcImage->imageData + srcImage->width * y ;
return (void**)rows;
}
}
void bhReleaseRowsA(void** rows)
{
free(rows);
}
void bhRemoveBorderA( BhImageA* srcImage ,BhImageA* dstImage,BhSizeA offset)
{
bhSetImageROIA(srcImage,bhRectA(offset.width,offset.height,dstImage->width,dstImage->height));
bhCopyImageA(srcImage,dstImage);
bhResetImageROIA(srcImage);
}
void bhSmoothExA(BhImageA* srcImage,BhImageA* dstImage,BhSizeA size,float sigma1,float sigma2,BhImageA* maskImage)
{
BhImageA* kernel = bhGet2DGaussianKernelA(size.width,sigma1);
BhSizeA border = bhSizeA(size.width /2,size.height /2);
BhImageA* srcImg = bhCreateFillBorderA(srcImage,border,BH_BORDER_CONSTANT,0);
BhImageA* dstImg = bhCreateBorderImageA(dstImage,border);
BhImageA* maskImg = bhCreateFillBorderA(maskImage,border,BH_BORDER_CONSTANT,0);
byte** maskRows = (byte**)bhCreateRowsA(maskImg);
BhRectA curBlock = bhRectA(0,0,size.width,size.height);
BhImageA* maskBlockImage = bhCreateImageHeaderA(size,8,maskImg->imageData);
maskBlockImage->widthStep = maskImage->widthStep;
if (srcImage->depth == BH_DEPTH_8)
{
float** kernelRows = (float**) bhCreateRowsA(kernel);
byte** dstRows = (byte**)bhCreateRowsA(dstImg);
byte** srcRows = (byte**)bhCreateRowsA(srcImg);
for (int i = border.height; i < srcImage->height - border.height; i++)
for (int j = border.width; j < srcImage->width - border.width; j++)
if (maskRows[i][j] != 0 )
{
curBlock.x = j - border.width;
curBlock.y = i - border.height;
maskBlockImage->imageData = &maskRows[i][j];
float kernelSum = bhSumA(kernel,maskBlockImage);
double sum = 0;
for (int n = -border.height; n <= border.height; n++)
for (int m = -border.width; m <= border.width; m++)
if (maskRows[n+i][m+j] == 255)
{
sum += ((double)srcRows[n+ i][m+j] * kernelRows[n+border.height][m+border.width] );
}
int k= 1;
if (kernelSum < 0.01)
MessageBoxA(0,"SS","",0);
else dstRows[i][j] = (byte)(sum /kernelSum );
}
bhReleaseRowsA((void**)kernelRows);
bhReleaseRowsA((void**)dstRows);
bhReleaseRowsA((void**)srcRows);
bhReleaseRowsA((void**)maskRows);
}
else if (srcImage->depth == BH_DEPTH_32_F)
{
float** kernelRows = (float**) bhCreateRowsA(kernel);
float** dstRows = (float**)bhCreateRowsA(dstImg);
float** srcRows = (float**)bhCreateRowsA(srcImg);
for (int i = border.height; i < srcImage->height - border.height; i++)
for (int j = border.width; j < srcImage->width - border.width; j++)
if (maskRows[i][j] != 0 )
{
curBlock.x = j - border.width;
curBlock.y = i - border.height;
maskBlockImage->imageData = &maskRows[i][j];
float kernelSum = bhSumA(kernel,maskBlockImage);
double sum = 0;
for (int n = -border.height; n <= border.height; n++)
for (int m = -border.width; m <= border.width; m++)
if (maskRows[n+i][m+j] == 255)
sum += ((double)srcRows[n+ i][m+j] * kernelRows[n+border.height][m+border.width] );
dstRows[i][j] = (sum /kernelSum );
}
bhReleaseRowsA((void**)kernelRows);
bhReleaseRowsA((void**)dstRows);
bhReleaseRowsA((void**)srcRows);
bhReleaseRowsA((void**)maskRows);
}
bhRemoveBorderA(dstImg,dstImage,border);
bhReleaseImageA(dstImg);
bhReleaseImageA(kernel);
bhReleaseImageA(srcImg);
bhReleaseImageA(maskImg);
bhReleaseImageHeaderA(maskBlockImage);
}