CPU DCT (OpenMP)

#include "dct_cpu.h"

#include <math.h>

#include <omp.h>

#ifndef M_PI

#define M_PI 3.14159265358979323846

#endif

#define M_PI_F ((float)M_PI)

// DCT normalization factors

static inline float alpha(int u) {

return (u == 0) ? (1.0f / sqrtf(2.0f)) : 1.0f;

// Forward 8x8 DCT-II

void forwardDCT8x8_CPU(const float* input, float* output) {

const int N = 8;

const float norm = 0.25f;

for (int v = 0; v < N; v++) {

for (int u = 0; u < N; u++) {

float sum = 0.0f;

 for (int y = 0; y < N; y++) {

for (int x = 0; x < N; x++) {

float px = input[y * N + x];

float cu = cosf((2.0f * x + 1.0f) * u * M_PI / (2.0f * N));

float cv = cosf((2.0f * y + 1.0f) * v * M_PI / (2.0f * N));

sum += px * cu * cv;

output[v * N + u] = norm * alpha(u) * alpha(v) * sum;

// Inverse 8x8 DCT-III

void inverseDCT8x8_CPU(const float* input, float* output) {

const int N = 8;

const float norm = 0.25f;

for (int y = 0; y < N; y++) {

for (int x = 0; x < N; x++) {

float sum = 0.0f;

for (int v = 0; v < N; v++) {

for (int u = 0; u < N; u++) {

float coeff = input[v * N + u];

 float cu = cosf((2.0f * x + 1.0f) * u * M_PI_F / (2.0f * N));

float cv = cosf((2.0f * y + 1.0f) * v * M_PI_F / (2.0f * N));

sum += alpha(u) * alpha(v) * coeff * cu * cv;

output[y * N + x] = norm * sum;

// Process entire image with block DCT

void blockDCT_CPU(const float* image, float* dctCoeffs,

int width, int height) {

const int BLOCK_SIZE = 8;

int numBlocksX = width / BLOCK_SIZE;

int numBlocksY = height / BLOCK_SIZE;

int by, bx;

#pragma omp parallel for private(by, bx) collapse(2)

for (by = 0; by < numBlocksY; by++) {

for (bx = 0; bx < numBlocksX; bx++) {

float block[64];

float dctBlock[64];

int y, x;
 // Extract 8x8 block

for (y = 0; y < BLOCK_SIZE; y++) {

for (x = 0; x < BLOCK_SIZE; x++) {

int imgY = by * BLOCK_SIZE + y;

int imgX = bx * BLOCK_SIZE + x;

block[y * BLOCK_SIZE + x] = image[imgY * width + imgX];

// Apply DCT

forwardDCT8x8_CPU(block, dctBlock);

// Store DCT block

for (y = 0; y < BLOCK_SIZE; y++) {

for (x = 0; x < BLOCK_SIZE; x++) {

int imgY = by * BLOCK_SIZE + y;

int imgX = bx * BLOCK_SIZE + x;

dctCoeffs[imgY * width + imgX] = dctBlock[y * BLOCK_SIZE + x];

// Process entire image with inverse block DCT

void inverseBlockDCT_CPU(const float* dctCoeffs, float* image,

int width, int height) {

const int BLOCK_SIZE = 8;

int numBlocksX = width / BLOCK_SIZE;

int numBlocksY = height / BLOCK_SIZE;

int by, bx;

#pragma omp parallel for private(by, bx) collapse(2)

for (by = 0; by < numBlocksY; by++) {

for (bx = 0; bx < numBlocksX; bx++) {

float block[64];

float dctBlock[64];

int y, x;

// Extract 8x8 DCT block

for (y = 0; y < BLOCK_SIZE; y++) {

for (x = 0; x < BLOCK_SIZE; x++) {

int imgY = by * BLOCK_SIZE + y;

int imgX = bx * BLOCK_SIZE + x;

dctBlock[y * BLOCK_SIZE + x] = dctCoeffs[imgY * width + imgX];

// Apply Inverse DCT

inverseDCT8x8_CPU(dctBlock, block);
 // Store image block

for (y = 0; y < BLOCK_SIZE; y++) {

for (x = 0; x < BLOCK_SIZE; x++) {

int imgY = by * BLOCK_SIZE + y;

int imgX = bx * BLOCK_SIZE + x;

image[imgY * width + imgX] = block[y * BLOCK_SIZE + x];

}
GPU DCT (CUDA)

#include "dct_cuda.cuh"

#include "utils.h"

#include <cuda_runtime.h>

#include <device_launch_parameters.h>

#define M_PI 3.14159265358979323846f

#define BLOCK_SIZE 8

// Device function for DCT normalization

__device__ inline float d_alpha(int u) {

return (u == 0) ? (1.0f / sqrtf(2.0f)) : 1.0f;

// Forward DCT kernel - each thread block processes one 8x8 image block

__global__ void forwardDCT8x8_kernel(float* blocks, int numBlocks, int width) {

__shared__ float sharedBlock[BLOCK_SIZE][BLOCK_SIZE];

__shared__ float sharedDCT[BLOCK_SIZE][BLOCK_SIZE];

int blockIdx_1d = blockIdx.y * gridDim.x + blockIdx.x;

if (blockIdx_1d >= numBlocks) return;

int tx = threadIdx.x;

int ty = threadIdx.y;
// Calculate block position in image

int blocksPerRow = width / BLOCK_SIZE;

int blockRow = blockIdx_1d / blocksPerRow;

int blockCol = blockIdx_1d % blocksPerRow;

int imgStartY = blockRow * BLOCK_SIZE;

int imgStartX = blockCol * BLOCK_SIZE;

// Load 8x8 block into shared memory

int imgY = imgStartY + ty;

int imgX = imgStartX + tx;

sharedBlock[ty][tx] = blocks[imgY * width + imgX];

__syncthreads();

// Compute DCT - each thread computes one coefficient

float sum = 0.0f;

const float norm = 0.25f;

for (int y = 0; y < BLOCK_SIZE; y++) {

for (int x = 0; x < BLOCK_SIZE; x++) {

float px = sharedBlock[y][x];

float cu = cosf((2.0f * x + 1.0f) * tx * M_PI / (2.0f * BLOCK_SIZE));

float cv = cosf((2.0f * y + 1.0f) * ty * M_PI / (2.0f * BLOCK_SIZE));

sum += px * cu * cv;
 }

sharedDCT[ty][tx] = norm * d_alpha(tx) * d_alpha(ty) * sum;

__syncthreads();

// Write DCT coefficients back to global memory

blocks[imgY * width + imgX] = sharedDCT[ty][tx];

// Inverse DCT kernel - each thread block processes one 8x8 DCT block

__global__ void inverseDCT8x8_kernel(float* blocks, int numBlocks, int width) {

__shared__ float sharedDCT[BLOCK_SIZE][BLOCK_SIZE];

__shared__ float sharedBlock[BLOCK_SIZE][BLOCK_SIZE];

int blockIdx_1d = blockIdx.y * gridDim.x + blockIdx.x;

if (blockIdx_1d >= numBlocks) return;

int tx = threadIdx.x;

int ty = threadIdx.y;

// Calculate block position

int blocksPerRow = width / BLOCK_SIZE;

int blockRow = blockIdx_1d / blocksPerRow;

int blockCol = blockIdx_1d % blocksPerRow;

int imgStartY = blockRow * BLOCK_SIZE;

int imgStartX = blockCol * BLOCK_SIZE;

// Load DCT block into shared memory

int imgY = imgStartY + ty;

int imgX = imgStartX + tx;

sharedDCT[ty][tx] = blocks[imgY * width + imgX];

__syncthreads();

// Compute inverse DCT - each thread computes one pixel

float sum = 0.0f;

const float norm = 0.25f;

for (int v = 0; v < BLOCK_SIZE; v++) {

for (int u = 0; u < BLOCK_SIZE; u++) {

float coeff = sharedDCT[v][u];

float cu = cosf((2.0f * tx + 1.0f) * u * M_PI / (2.0f * BLOCK_SIZE));

float cv = cosf((2.0f * ty + 1.0f) * v * M_PI / (2.0f * BLOCK_SIZE));

sum += d_alpha(u) * d_alpha(v) * coeff * cu * cv;

sharedBlock[ty][tx] = norm * sum;

 __syncthreads();

// Write reconstructed block back to global memory

blocks[imgY * width + imgX] = sharedBlock[ty][tx];

// Host wrapper for forward DCT

void forwardDCT_CUDA(float* d_image, int width, int height) {

int numBlocksX = width / BLOCK_SIZE;

int numBlocksY = height / BLOCK_SIZE;

int totalBlocks = numBlocksX * numBlocksY;

// Each CUDA block processes one 8x8 image block

dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);

// Grid dimensions

int gridX = (totalBlocks > 65535) ? 65535 : totalBlocks;

int gridY = (totalBlocks + gridX - 1) / gridX;

dim3 numBlocks(gridX, gridY);

forwardDCT8x8_kernel<<<numBlocks, threadsPerBlock>>>(d_image, totalBlocks, width);

cudaCheckError(cudaGetLastError());

cudaCheckError(cudaDeviceSynchronize());

// Host wrapper for inverse DCT

void inverseDCT_CUDA(float* d_image, int width, int height) {

int numBlocksX = width / BLOCK_SIZE;

int numBlocksY = height / BLOCK_SIZE;

int totalBlocks = numBlocksX * numBlocksY;

dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);

int gridX = (totalBlocks > 65535) ? 65535 : totalBlocks;

int gridY = (totalBlocks + gridX - 1) / gridX;

dim3 numBlocks(gridX, gridY);

inverseDCT8x8_kernel<<<numBlocks, threadsPerBlock>>>(d_image, totalBlocks, width);

cudaCheckError(cudaGetLastError());

cudaCheckError(cudaDeviceSynchronize());

}
Image I/O

#define _CRT_SECURE_NO_WARNINGS

#include "image_io.h"

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <cuda_runtime.h>

#include "utils.h" // For cudaCheckError

#define STB_IMAGE_IMPLEMENTATION

#include "../include/stb_image.h"

#define STB_IMAGE_WRITE_IMPLEMENTATION

#include "../include/stb_image_write.h"

// Helper for clamping (defined as macro in utils.h)

Image* loadImage(const char* filename) {

int width, height, channels;

// Force 3 channels (RGB)

unsigned char* data = stbi_load(filename, &width, &height, &channels, 3);

if (!data) {

fprintf(stderr, "Error: Could not load image %s\n", filename);

 return NULL;

Image* img = (Image*)malloc(sizeof(Image));

if (!img) {

stbi_image_free(data);

return NULL;

img->width = width;

img->height = height;

img->channels = 3; // Forced to 3

img->data = data;

printf("Loaded image: %s (%dx%d, %d channels)\n",

filename, img->width, img->height, img->channels);

return img;

void freeImage(Image* img) {

if (img) {

if (img->data) {

stbi_image_free(img->data);

free(img);
 }

int saveImage(const char* filename, Image* img) {

int result = 0;

// Determine file format from extension

// Simple check: default to PNG if unknown or png

// strstr is case sensitive, but for simplicity we check lowercase commonly

// In C, no std::string::find

const char* ext = strrchr(filename, '.');

int use_jpg = 0;

int use_bmp = 0;

if (ext) {

if (strstr(ext, "jpg") || strstr(ext, "JPG") || strstr(ext, "jpeg")) use_jpg = 1;

else if (strstr(ext, "bmp") || strstr(ext, "BMP")) use_bmp = 1;

if (use_jpg) {

result = stbi_write_jpg(filename, img->width, img->height,

img->channels, img->data, 95); // Quality 95

} else if (use_bmp) {

result = stbi_write_bmp(filename, img->width, img->height,

img->channels, img->data);
 } else {

// Default PNG

result = stbi_write_png(filename, img->width, img->height,

img->channels, img->data, img->width * img->channels);

if (result) {

printf("Saved image: %s\n", filename);

return 1;

} else {

fprintf(stderr, "Error: Could not save image %s\n", filename);

return 0;

void rgb2ycbcr(Image* img) {

if (img->channels != 3) {

fprintf(stderr, "Warning: rgb2ycbcr only works on 3-channel images\n");

return;

int totalPixels = img->width * img->height;

for (int i = 0; i < totalPixels; i++) {

int idx = i * 3;

float r = img->data[idx + 0];

 float g = img->data[idx + 1];

float b = img->data[idx + 2];

// ITU-R BT.601 conversion

float y = 0.299f * r + 0.587f * g + 0.114f * b;

float cb = -0.168736f * r - 0.331264f * g + 0.5f * b + 128.0f;

float cr = 0.5f * r - 0.418688f * g - 0.081312f * b + 128.0f;

img->data[idx + 0] = (unsigned char)clamp(y, 0.0f, 255.0f);

img->data[idx + 1] = (unsigned char)clamp(cb, 0.0f, 255.0f);

img->data[idx + 2] = (unsigned char)clamp(cr, 0.0f, 255.0f);

printf("Converting to YCbCr color space...\n");

void ycbcr2rgb(Image* img) {

if (img->channels != 3) {

fprintf(stderr, "Warning: ycbcr2rgb only works on 3-channel images\n");

return;

int totalPixels = img->width * img->height;

for (int i = 0; i < totalPixels; i++) {

int idx = i * 3;

float y = img->data[idx + 0];

 float cb = img->data[idx + 1] - 128.0f;

float cr = img->data[idx + 2] - 128.0f;

// ITU-R BT.601 inverse conversion

float r = y + 1.402f * cr;

float g = y - 0.344136f * cb - 0.714136f * cr;

float b = y + 1.772f * cb;

img->data[idx + 0] = clamp(r, 0.0f, 255.0f);

img->data[idx + 1] = clamp(g, 0.0f, 255.0f);

img->data[idx + 2] = clamp(b, 0.0f, 255.0f);

printf("Converting back to RGB...\n");

// Bilinear interpolation resize

Image* resizeImage(Image* src, int targetWidth, int targetHeight) {

Image* dst = (Image*)malloc(sizeof(Image));

if (!dst) return NULL;

dst->width = targetWidth;

dst->height = targetHeight;

dst->channels = src->channels;

dst->data = (unsigned char*)malloc(targetWidth * targetHeight * src->channels);

if (!dst->data) {
 free(dst);

return NULL;

float x_ratio = ((float)(src->width - 1)) / targetWidth;

float y_ratio = ((float)(src->height - 1)) / targetHeight;

for (int i = 0; i < targetHeight; i++) {

for (int j = 0; j < targetWidth; j++) {

int x = (int)(x_ratio * j);

int y = (int)(y_ratio * i);

float x_diff = (x_ratio * j) - x;

float y_diff = (y_ratio * i) - y;

int index = (y * src->width + x) * src->channels;

int b_index = (i * targetWidth + j) * src->channels;

// For each channel

for (int k = 0; k < src->channels; k++) {

// Get pixels at (x,y), (x+1,y), (x,y+1), (x+1,y+1)

float A = (float)src->data[index + k];

float B = (float)src->data[index + src->channels + k];

float C = (float)src->data[index + src->width * src->channels + k];

float D = (float)src->data[index + src->width * src->channels + src->channels + k];

// Bilinear interpolation formula

 float val = A * (1 - x_diff) * (1 - y_diff) +

B * (x_diff) * (1 - y_diff) +

C * (y_diff) * (1 - x_diff) +

D * (x_diff * y_diff);

dst->data[b_index + k] = clamp(val, 0.0f, 255.0f);

printf("Resized image from %dx%d to %dx%d\n", src->width, src->height, targetWidth,

targetHeight);

return dst;

float* allocateDeviceImage(int width, int height, int channels) {

int size = width * height * channels;

float* d_img;

cudaCheckError(cudaMalloc((void**)&d_img, size * sizeof(float)));

return d_img;

void copyImageToDevice(float* d_img, const unsigned char* h_img, int size) {

float* h_float = (float*)malloc(size * sizeof(float));

 if (!h_float) {

fprintf(stderr, "Error: Out of memory in copyImageToDevice\n");

return;

for (int i = 0; i < size; i++) {

h_float[i] = (float)h_img[i];

cudaCheckError(cudaMemcpy(d_img, h_float, size * sizeof(float),

cudaMemcpyHostToDevice));

free(h_float);

void copyImageToHost(unsigned char* h_img, const float* d_img, int size) {

float* h_float = (float*)malloc(size * sizeof(float));

if (!h_float) {

fprintf(stderr, "Error: Out of memory in copyImageToHost\n");

return;

cudaCheckError(cudaMemcpy(h_float, d_img, size * sizeof(float),

cudaMemcpyDeviceToHost));

for (int i = 0; i < size; i++) {

h_img[i] = clamp(h_float[i], 0.0f, 255.0f);

 }

free(h_float);

void freeDeviceImage(float* d_img) {

if (d_img) {

cudaCheckError(cudaFree(d_img));

}
CLI interface

#define _CRT_SECURE_NO_WARNINGS

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "utils.h"

#include "image_io.h"

#include "dct_cpu.h"

#include "dct_cuda.cuh"

#include "watermark_cox.h"

#include "watermark_cox_cuda.cuh"

// Helper to check command line flags

int check_arg(const char* arg, const char* long_opt, const char* short_opt) {

if (strcmp(arg, long_opt) == 0) return 1;

if (short_opt && strcmp(arg, short_opt) == 0) return 1;

return 0;

void printUsage() {

printf("\n=== Cox Spread Spectrum Watermarking System ===\n");

printf("\nUsage:\n");

printf(" ImageWatermarking --keygen -k <keyfile> --length <k> [--alpha <alpha>] [--seed

<seed>]\n");

printf(" ImageWatermarking --embed -i <input> -o <output> -k <keyfile> [--cpu|--gpu]\n");

 printf(" ImageWatermarking --extract -i <watermarked> -k <keyfile> -o <original> [--cpu|--
gpu]\n");

printf(" ImageWatermarking --benchmark -i <input> -k <keyfile>\n");

printf(" ImageWatermarking --gpuinfo\n");

printf("\nOptions:\n");

printf(" --keygen Generate watermark key\n");

printf(" --embed Embed watermark into image\n");

printf(" --extract Extract/detect watermark from image\n");

printf(" --benchmark Benchmark CPU vs GPU performance\n");

printf(" --gpuinfo Display GPU information\n");

printf(" -i, --input Input image file\n");

printf(" -o, --output Output image file\n");

printf(" -k, --key Watermark key file\n");

printf(" --length Number of coefficients (default: 1000)\n");

printf(" --alpha Embedding strength (default: 0.1)\n");

printf(" --seed Random seed (default: 12345)\n");

printf(" --cpu Use CPU implementation\n");

printf(" --gpu Use GPU implementation (default)\n");

printf("\nExamples:\n");

printf(" ImageWatermarking --keygen -k [Link] --length 1000 --alpha 0.1\n");

printf(" ImageWatermarking --embed -i [Link] -o [Link] -k [Link]\n");

printf(" ImageWatermarking --extract -i [Link] -k [Link] -o [Link]\n");

printf("\n");

int main(int argc, char** argv) {

 if (argc < 2) {

printUsage();

return 1;

const char* mode = NULL;

const char* inputFile = NULL;

const char* outputFile = NULL;

const char* keyFile = NULL;

const char* originalFile = NULL; // Not strictly needed if input/output logic handles it but
helpful

int k = 1000;

float alpha = 0.1f;

unsigned int seed = 12345;

int useGPU = 1;

// Parse arguments

int i;

for (i = 1; i < argc; i++) {

if (check_arg(argv[i], "--keygen", NULL)) mode = "keygen";

else if (check_arg(argv[i], "--embed", NULL)) mode = "embed";

else if (check_arg(argv[i], "--extract", NULL)) mode = "extract";

else if (check_arg(argv[i], "--benchmark", NULL)) mode = "benchmark";

else if (check_arg(argv[i], "--gpuinfo", NULL)) mode = "gpuinfo";

else if (check_arg(argv[i], "--input", "-i")) inputFile = argv[++i];

 else if (check_arg(argv[i], "--output", "-o")) outputFile = argv[++i];

else if (check_arg(argv[i], "--key", "-k")) keyFile = argv[++i];

else if (check_arg(argv[i], "--length", NULL)) k = atoi(argv[++i]);

else if (check_arg(argv[i], "--alpha", NULL)) alpha = (float)atof(argv[++i]);

else if (check_arg(argv[i], "--seed", NULL)) seed = (unsigned int)strtoul(argv[++i], NULL,

10);

else if (check_arg(argv[i], "--cpu", NULL)) useGPU = 0;

else if (check_arg(argv[i], "--gpu", NULL)) useGPU = 1;

if (!mode) {

printUsage();

return 1;

// Execute mode

if (strcmp(mode, "gpuinfo") == 0) {

printGPUInfo();

return 0;

if (strcmp(mode, "keygen") == 0) {

if (!keyFile) {

fprintf(stderr, "Error: Key file not specified\n");

return 1;

}
 printf("Generating watermark key...\n");

printf(" Length (k): %d\n", k);

printf(" Alpha: %.2f\n", alpha);

printf(" Seed: %u\n", seed);

WatermarkKey* key = generateKey(k, seed, alpha);

if (key) {

saveKey(keyFile, key);

freeKey(key);

printf("Key generation complete!\n");

return 0;

if (strcmp(mode, "embed") == 0) {

if (!inputFile || !outputFile || !keyFile) {

fprintf(stderr, "Error: Input, output, or key file missing\n");

return 1;

printf("\n=== Watermark Embedding ===\n");

printf("Mode: %s\n", useGPU ? "GPU (CUDA)" : "CPU");

WatermarkKey* key = loadKey(keyFile);

if (!key) return 1;
Image* img = loadImage(inputFile);

if (!img) {

freeKey(key);

return 1;

rgb2ycbcr(img);

int width = img->width;

int height = img->height;

int channels = img->channels; // 3

int size = width * height;

// Prepare Y channel for DCT

float* yChannel = (float*)malloc(size * sizeof(float));

if (!yChannel) return 1;

int i;

for (i = 0; i < size; i++) {

yChannel[i] = (float)img->data[i * 3]; // Y component

float* dctCoeffs = (float*)malloc(size * sizeof(float));

 Timer timer;

initTimer(&timer);

startTimer(&timer);

if (useGPU) {

// Alloc GPU mem

float* d_img = allocateDeviceImage(width, height, 1);

float* d_dct = allocateDeviceImage(width, height, 1); // Unused logic kept for

structure

// Copy Y to GPU

cudaCheckError(cudaMemcpy(d_img, yChannel, size * sizeof(float),

cudaMemcpyHostToDevice));

// DCT

forwardDCT_CUDA(d_img, width, height);

// Calculate indices if needed (Host side)

if (key->indices[0] == -1) {

float* h_dct = (float*)malloc(size * sizeof(float));

if (h_dct) {

cudaCheckError(cudaMemcpy(h_dct, d_img, size * sizeof(float),

cudaMemcpyDeviceToHost));

selectCoefficients(h_dct, key->indices, key->k, width, height);

free(h_dct);

}
 }

// Embed

embedWatermark_Cox_CUDA(d_img, key);

// IDCT

inverseDCT_CUDA(d_img, width, height);

// Copy back

cudaCheckError(cudaMemcpy(yChannel, d_img, size * sizeof(float),

cudaMemcpyDeviceToHost));

freeDeviceImage(d_img);

freeDeviceImage(d_dct);

} else {

// CPU Mode

// DCT

blockDCT_CPU(yChannel, dctCoeffs, width, height);

// Embed

embedWatermark_Cox_CPU(dctCoeffs, key, width, height);

// IDCT

inverseBlockDCT_CPU(dctCoeffs, yChannel, width, height);

}
{

double elapsed = stopTimer(&timer);

if (useGPU) printf("GPU processing time: %.3f ms\n", elapsed);

else printf("CPU processing time: %.3f ms\n", elapsed);

// Put Y back into image

int i;

for (i = 0; i < size; i++) {

img->data[i * 3] = (unsigned char)clamp(yChannel[i], 0.0f, 255.0f);

ycbcr2rgb(img);

if (saveImage(outputFile, img)) {

// Image saved

printf("Watermark embedding complete!\n");

free(yChannel);

free(dctCoeffs); // Only used in CPU path but safe to free if null? No, malloced.

freeImage(img);

freeKey(key);
 return 0;

if (strcmp(mode, "extract") == 0) {

if (!inputFile || !outputFile || !keyFile) {

if (!outputFile) {

fprintf(stderr, "Error: Original reference image required (-o)\n");

return 1;

originalFile = outputFile; // Remap

const char* watermarkedFile = inputFile;

printf("\n=== Watermark Extraction ===\n");

printf("Mode: %s\n", useGPU ? "GPU (CUDA)" : "CPU");

WatermarkKey* key = loadKey(keyFile);

if (!key) return 1;

Image* imgWM = loadImage(watermarkedFile);

Image* imgOrig = loadImage(originalFile);

if (!imgWM || !imgOrig) return 1;

rgb2ycbcr(imgWM);

rgb2ycbcr(imgOrig);

// Validation check

if (imgWM->width != imgOrig->width || imgWM->height != imgOrig->height) {

printf("Warning: Image dimensions mismatch!\n");

printf(" Watermarked: %dx%d\n", imgWM->width, imgWM->height);

printf(" Original: %dx%d\n", imgOrig->width, imgOrig->height);

printf("Resizing watermarked image to match original...\n");

Image* resizedWM = resizeImage(imgWM, imgOrig->width, imgOrig->height);

if (!resizedWM) return 1;

freeImage(imgWM);

imgWM = resizedWM;

int size = imgOrig->width * imgOrig->height;

float* yWM = (float*)malloc(size * sizeof(float));

float* yOrig = (float*)malloc(size * sizeof(float));

int i;

for(i=0; i<size; i++) {

yWM[i] = (float)imgWM->data[i*3];

yOrig[i] = (float)imgOrig->data[i*3];
 }

float correlation = 0.0f;

Timer timer;

initTimer(&timer);

startTimer(&timer);

if (useGPU) {

float* d_wm = allocateDeviceImage(imgOrig->width, imgOrig->height, 1);

float* d_orig = allocateDeviceImage(imgOrig->width, imgOrig->height, 1);

cudaCheckError(cudaMemcpy(d_wm, yWM, size*sizeof(float),

cudaMemcpyHostToDevice));

cudaCheckError(cudaMemcpy(d_orig, yOrig, size*sizeof(float),

cudaMemcpyHostToDevice));

forwardDCT_CUDA(d_wm, imgOrig->width, imgOrig->height);

forwardDCT_CUDA(d_orig, imgOrig->width, imgOrig->height);

// Calculate indices from ORIGINAL image if needed

if (key->indices[0] == -1) {

float* h_dct = (float*)malloc(size * sizeof(float));

if (h_dct) {

cudaCheckError(cudaMemcpy(h_dct, d_orig, size * sizeof(float),

cudaMemcpyDeviceToHost));

selectCoefficients(h_dct, key->indices, key->k, imgOrig->width, imgOrig->height);

 free(h_dct);

correlation = extractWatermark_Cox_CUDA(d_orig, d_wm, key);

freeDeviceImage(d_wm);

freeDeviceImage(d_orig);

} else {

float* dctWM = (float*)malloc(size*sizeof(float));

float* dctOrig = (float*)malloc(size*sizeof(float));

blockDCT_CPU(yWM, dctWM, imgOrig->width, imgOrig->height);

blockDCT_CPU(yOrig, dctOrig, imgOrig->width, imgOrig->height);

correlation = extractWatermark_Cox_CPU(dctOrig, dctWM, key, imgOrig->width,

imgOrig->height);

free(dctWM);

free(dctOrig);

double elapsed = stopTimer(&timer);

printf("Watermark correlation: %.4f\n", correlation);

if (correlation > 6.0f) printf("CHECKMARK Watermark DETECTED (strong

correlation)\n");
 else if (correlation > 3.0f) printf("Warning: Weak watermark detection\n");

else printf("CROSS Watermark NOT detected or destroyed\n");

if (useGPU) printf("GPU processing time: %.3f ms\n", elapsed);

else printf("CPU processing time: %.3f ms\n", elapsed);

printf("Watermark extraction complete!\n");

printf("Final correlation: %.4f\n", correlation);

free(yWM);

free(yOrig);

freeImage(imgWM);

freeImage(imgOrig);

freeKey(key);

return 0;

if (strcmp(mode, "benchmark") == 0) {

if (!inputFile || !keyFile) {

fprintf(stderr, "Error: Input image and key file required for benchmark\n");

return 1;

printf("\n=== Performance Benchmark Suite ===\n");

printf("Testing: %s\n\n", inputFile);

WatermarkKey* key = loadKey(keyFile);

if (!key) return 1;

Image* img = loadImage(inputFile);

if (!img) {

freeKey(key);

return 1;

rgb2ycbcr(img);

int width = img->width;

int height = img->height;

int size = width * height;

float* yChannel = (float*)malloc(size * sizeof(float));

float* dctCoeffs = (float*)malloc(size * sizeof(float));

if (!yChannel || !dctCoeffs) return 1;

int i;

for (i = 0; i < size; i++) {

yChannel[i] = (float)img->data[i * 3];

printf("Image Size: %dx%d (%d pixels)\n", width, height, size);

printf("Watermark Coefficients: %d\n", key->k);

printf("Alpha: %.2f\n\n", key->alpha);

// CPU Benchmark

printf("--- CPU Performance ---\n");

Timer cpuTimer;

initTimer(&cpuTimer);

startTimer(&cpuTimer);

blockDCT_CPU(yChannel, dctCoeffs, width, height);

embedWatermark_Cox_CPU(dctCoeffs, key, width, height);

inverseBlockDCT_CPU(dctCoeffs, yChannel, width, height);

double cpuTime = stopTimer(&cpuTimer);

printf("Total CPU Time: %.3f ms\n\n", cpuTime);

// Reset data

for (i = 0; i < size; i++) {

yChannel[i] = (float)img->data[i * 3];

// GPU Benchmark

printf("--- GPU Performance ---\n");

float* d_img = allocateDeviceImage(width, height, 1);

Timer gpuTimer;

initTimer(&gpuTimer);
 startTimer(&gpuTimer);

cudaCheckError(cudaMemcpy(d_img, yChannel, size * sizeof(float),

cudaMemcpyHostToDevice));

forwardDCT_CUDA(d_img, width, height);

if (key->indices[0] == -1) {

float* h_dct = (float*)malloc(size * sizeof(float));

if (h_dct) {

cudaCheckError(cudaMemcpy(h_dct, d_img, size * sizeof(float),

cudaMemcpyDeviceToHost));

selectCoefficients(h_dct, key->indices, key->k, width, height);

free(h_dct);

embedWatermark_Cox_CUDA(d_img, key);

inverseDCT_CUDA(d_img, width, height);

cudaCheckError(cudaMemcpy(yChannel, d_img, size * sizeof(float),

cudaMemcpyDeviceToHost));

double gpuTime = stopTimer(&gpuTimer);

printf("Total GPU Time: %.3f ms\n\n", gpuTime);

freeDeviceImage(d_img);

// Results Summary
 double speedup = cpuTime / gpuTime;

printf("=== Benchmark Results ===\n");

printf("CPU Time: %.3f ms\n", cpuTime);

printf("GPU Time: %.3f ms\n", gpuTime);

printf("Speedup: %.2fx\n", speedup);

printf("Efficiency: %.1f%%\n\n", (speedup / 1.0) * 100.0);

// Export to CSV

FILE* csvFile = fopen("benchmark_results.csv", "a");

if (csvFile) {

fseek(csvFile, 0, SEEK_END);

if (ftell(csvFile) == 0) {

fprintf(csvFile,
"ImageSize,Width,Height,Pixels,Coefficients,CPU_ms,GPU_ms,Speedup\n");

fprintf(csvFile, "%dx%d,%d,%d,%d,%d,%.3f,%.3f,%.2f\n",

width, height, width, height, size, key->k, cpuTime, gpuTime, speedup);

fclose(csvFile);

printf("Results appended to benchmark_results.csv\n");

free(yChannel);

free(dctCoeffs);

freeImage(img);

freeKey(key);

return 0;
 }

return 0;

}
GPU watermarking

#include "watermark_cox_cuda.cuh"

#include "utils.h"

#include <cuda_runtime.h>

#include <device_launch_parameters.h>

#include <stdio.h>

#include <math.h>

// CUDA kernel for watermark embedding

__global__ void embedWatermark_Cox_kernel(float* dctCoeffs,

const int* indices,

const float* watermark,

int k, float alpha) {

int i = blockIdx.x * blockDim.x + threadIdx.x;

if (i < k) {

int idx = indices[i];

float V = dctCoeffs[idx];

float W = watermark[i];

// Cox formula: V' = V * (1 + alpha * W)

dctCoeffs[idx] = V * (1.0f + alpha * W);

}
// CUDA kernel for watermark extraction using atomic adds

__global__ void extractWatermark_Cox_kernel(const float* dctCoeffsOrig,

const float* dctCoeffsTest,

const int* indices,

const float* watermark,

int k,

float* partialSums) {

int i = blockIdx.x * blockDim.x + threadIdx.x;

if (i < k) {

int idx = indices[i];

float V = dctCoeffsOrig[idx];

float V_star = dctCoeffsTest[idx];

float W = watermark[i];

float diff = 0.0f;

if (fabsf(V) > 1e-6f) {

diff = (V_star - V) / V;

float num = diff * W;

float den = diff * diff;

atomicAdd(&partialSums[0], num); // Accumulate numerator

atomicAdd(&partialSums[1], den); // Accumulate denominator

}
}

// Host wrapper for embedding

void embedWatermark_Cox_CUDA(float* d_dctCoeffs, WatermarkKey* key) {

// Copy key data to device

int* d_indices;

float* d_watermark;

cudaCheckError(cudaMalloc((void**)&d_indices, key->k * sizeof(int)));

cudaCheckError(cudaMalloc((void**)&d_watermark, key->k * sizeof(float)));

cudaCheckError(cudaMemcpy(d_indices, key->indices, key->k * sizeof(int),

cudaMemcpyHostToDevice));

cudaCheckError(cudaMemcpy(d_watermark, key->watermark, key->k * sizeof(float),

cudaMemcpyHostToDevice));

// Launch kernel

int blockSize = 256;

int numBlocks = (key->k + blockSize - 1) / blockSize;

embedWatermark_Cox_kernel<<<numBlocks, blockSize>>>(

d_dctCoeffs, d_indices, d_watermark, key->k, key->alpha

);

cudaCheckError(cudaGetLastError());

cudaCheckError(cudaDeviceSynchronize());
 // Cleanup

cudaCheckError(cudaFree(d_indices));

cudaCheckError(cudaFree(d_watermark));

printf("Watermark embedded using CUDA (alpha=%.2f)\n", key->alpha);

// Host wrapper for extraction

float extractWatermark_Cox_CUDA(const float* d_dctCoeffsOrig,

const float* d_dctCoeffsTest,

WatermarkKey* key) {

// Copy key data to device

int* d_indices;

float* d_watermark;

float* d_sums;

float h_sums[2] = {0.0f, 0.0f};

cudaCheckError(cudaMalloc((void**)&d_indices, key->k * sizeof(int)));

cudaCheckError(cudaMalloc((void**)&d_watermark, key->k * sizeof(float)));

cudaCheckError(cudaMalloc((void**)&d_sums, 2 * sizeof(float)));

cudaCheckError(cudaMemcpy(d_indices, key->indices, key->k * sizeof(int),

cudaMemcpyHostToDevice));

cudaCheckError(cudaMemcpy(d_watermark, key->watermark, key->k * sizeof(float),

cudaMemcpyHostToDevice));
 cudaCheckError(cudaMemcpy(d_sums, h_sums, 2 * sizeof(float),
cudaMemcpyHostToDevice));

// Launch kernel

int blockSize = 256;

int numBlocks = (key->k + blockSize - 1) / blockSize;

extractWatermark_Cox_kernel<<<numBlocks, blockSize>>>(

d_dctCoeffsOrig, d_dctCoeffsTest, d_indices, d_watermark, key->k, d_sums

);

cudaCheckError(cudaGetLastError());

cudaCheckError(cudaDeviceSynchronize());

// Copy results back

cudaCheckError(cudaMemcpy(h_sums, d_sums, 2 * sizeof(float),

cudaMemcpyDeviceToHost));

float correlation = 0.0f;

if (h_sums[1] > 1e-9f) {

correlation = h_sums[0] / sqrtf(h_sums[1]);

printf("Watermark correlation (CUDA): %.4f\n", correlation);

if (correlation > 6.0f) {

printf("CHECKMARK Watermark DETECTED (strong correlation)\n");

 } else if (correlation > 3.0f) {

printf("WARNING Watermark possibly detected (weak correlation)\n");

} else {

printf("CROSS Watermark NOT detected or destroyed\n");

// Cleanup

cudaCheckError(cudaFree(d_indices));

cudaCheckError(cudaFree(d_watermark));

cudaCheckError(cudaFree(d_sums));

return correlation;

}
watermark_cox.c

/**

* @file watermark_cox.c

* @brief CPU implementation of Cox spread spectrum watermarking

* Implements the Cox et al. (1997) algorithm for DCT-based watermarking.

* Parallelized using OpenMP for multi-core CPU acceleration.

* Algorithm:

* - Embedding: V'[i] = V[i] * (1 + alpha * W[i])

* - Detection: Normalized correlation coefficient

*/

#include "watermark_cox.h"

#include <math.h>

#include <stdlib.h>

#include <stdio.h>

// Structure for sorting coefficients by magnitude

typedef struct {

int index;

float magnitude;

} CoeffMag;
// Comparator for qsort (Descending order of magnitude, ascending order of index for
stability)

int compareCoeffMag(const void* a, const void* b) {

CoeffMag* cmA = (CoeffMag*)a;

CoeffMag* cmB = (CoeffMag*)b;

if (cmA->magnitude < cmB->magnitude) return 1;

if (cmA->magnitude > cmB->magnitude) return -1;

// Tie-breaker: sort by index (ascending) to ensure stability

if (cmA->index < cmB->index) return -1;

if (cmA->index > cmB->index) return 1;

return 0;

// Select k largest magnitude coefficients (excluding DC and very high frequencies)

void selectCoefficients(const float* dctCoeffs, int* indices, int k,

int width, int height) {

const int BLOCK_SIZE = 8;

int numBlocksX = width / BLOCK_SIZE;

int numBlocksY = height / BLOCK_SIZE;

// Estimate max possible candidates to allocate memory

// Each block has 64 coeffs. We skip some. Safe upper bound is total pixels.

int maxCandidates = width * height;

CoeffMag* coeffs = (CoeffMag*)malloc(maxCandidates * sizeof(CoeffMag));

if (!coeffs) {

fprintf(stderr, "Error: Out of memory in selectCoefficients\n");

return;

int count = 0;

int by, bx;

// Collect all mid-frequency coefficients

for (by = 0; by < numBlocksY; by++) {

for (bx = 0; bx < numBlocksX; bx++) {

int y, x;

for (y = 0; y < BLOCK_SIZE; y++) {

for (x = 0; x < BLOCK_SIZE; x++) {

int imgY, imgX, idx;

// Skip DC coefficient (0,0) and very high frequencies (7,7)

if ((x == 0 && y == 0) || (x >= 6 && y >= 6)) continue;

imgY = by * BLOCK_SIZE + y;

imgX = bx * BLOCK_SIZE + x;

idx = imgY * width + imgX;

coeffs[count].index = idx;

coeffs[count].magnitude = fabsf(dctCoeffs[idx]);

count++;
 }

// Sort entire array by magnitude (using qsort)

// Note: C++ partial_sort is faster but qsort is standard C

qsort(coeffs, count, sizeof(CoeffMag), compareCoeffMag);

// Select top k

int i;

for (i = 0; i < k && i < count; i++) {

indices[i] = coeffs[i].index;

free(coeffs);

printf("Selected %d largest DCT coefficients for watermarking\n", k);

// Embed watermark using Cox et al. (1997) formula: V'[i] = V[i] * (1 + alpha * W[i])

void embedWatermark_Cox_CPU(float* dctCoeffs, WatermarkKey* key,

int width, int height) {

int i;
 // If indices not set, select them now

if (key->indices[0] == -1) {

selectCoefficients(dctCoeffs, key->indices, key->k, width, height);

// Embed watermark

for (i = 0; i < key->k; i++) {

int idx = key->indices[i];

float V = dctCoeffs[idx];

float W = key->watermark[i];

// Cox formula: V' = V * (1 + alpha * W)

dctCoeffs[idx] = V * (1.0f + key->alpha * W);

printf("Watermark embedded using Cox algorithm (alpha=%.2f)\n", key->alpha);

// Extract and detect watermark using normalized correlation

float extractWatermark_Cox_CPU(const float* dctCoeffsOrig,

const float* dctCoeffsTest,

WatermarkKey* key,

int width, int height) {

float numerator = 0.0f;

 float denominator = 0.0f;

int i;

for (i = 0; i < key->k; i++) {

int idx = key->indices[i];

float V = dctCoeffsOrig[idx];

float V_star = dctCoeffsTest[idx];

float W = key->watermark[i];

float diff = 0.0f;

if (fabsf(V) > 1e-6f) {

diff = (V_star - V) / V;

// Normalized correlation: rho = sum(ext_W * ref_W) / sqrt(sum(ext_W^2))

numerator += diff * W;

denominator += diff * diff;

if (denominator < 1e-9f) return 0.0f;

return numerator / sqrtf(denominator);