图像处理 -- 仿射变换之Affine Transformation

仿射变换（Affine Transformation）

仿射变换是图像处理中的一种基本操作，通过线性变换和平移实现图像的几何变换。仿射变换包括旋转、缩放、平移、翻转、错切（shear）等操作。

1. 仿射变换的作用

• 旋转：将图像绕一个固定点（通常是图像中心）旋转一定角度。
• 缩放：对图像进行放大或缩小，改变图像的尺寸。
• 平移：将图像在平面上移动到指定位置。
• 翻转：沿特定轴将图像反转，比如水平或垂直翻转。
• 错切：将图像的某一个方向拉伸，使图像产生倾斜效果。

仿射变换在计算机视觉和图像处理中的应用非常广泛，比如在图像配准、目标检测、图像增强等领域中都可以看到仿射变换的身影。

2. 数学实现原理

仿射变换通过一个矩阵乘法和一个向量加法来实现。假设二维空间中的一个点的坐标为 $(x,y)$，经过仿射变换后，该点的新坐标为 $(x^{\prime },y^{\prime })$。则有如下数学表达式：

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}e\\ f\end{array}\right)$$

这里，$\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$ 是 2x2 仿射变换矩阵，$\left(\begin{array}{c}e\\ f\end{array}\right)$ 是平移向量。

具体仿射变换：

1. 平移：
   $$x^{\prime }=x+e,y^{\prime }=y+f$$
   其中，$e$ 和 $f$ 是平移量。

2. 旋转：
   $$x^{\prime }=x\cdot \cos (\theta )-y\cdot \sin (\theta ),y^{\prime }=x\cdot \sin (\theta )+y\cdot \cos (\theta )$$
   其中，$\theta$ 是旋转角度。

3. 缩放：
   $$x^{\prime }=x\cdot s_x,y^{\prime }=y\cdot s_y$$
   其中，$s_x$ 和 $s_y$ 是在 x 轴和 y 轴上的缩放比例。

4. 错切：
   $$x^{\prime }=x+{\lambda }_y\cdot y,y^{\prime }=y+{\lambda }_x\cdot x$$
   其中，${\lambda }_x$ 和 ${\lambda }_y$ 是在 x 轴和 y 轴上的错切系数。

3. 注意事项

1. 矩阵可逆性：为了保证仿射变换能被逆变换，仿射变换矩阵 $\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$ 必须是非奇异的，即其行列式不为 0。这确保了变换后的图像可以通过逆变换恢复原状。

2. 边界处理：仿射变换可能导致图像的一部分移出图像边界，或者图像空白区域增加。因此，通常需要对图像进行裁剪或填充处理。

3. 插值方法：在仿射变换过程中，目标图像中的像素可能对应源图像中非整数坐标的点，因此需要使用插值方法（如最近邻插值、双线性插值或双三次插值）来计算这些点的像素值。

4. 保持图像质量：频繁进行仿射变换可能导致图像质量下降，比如模糊、锯齿效应等。因此，在设计变换时需要考虑如何最小化质量损失。

5. 变换顺序：如果需要执行多个变换，变换的顺序会影响最终结果。例如，旋转后再平移和先平移后旋转的结果是不一样的。因此，需要谨慎设计变换的顺序以达到预期效果。

C代码示例

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define IMAGE_WIDTH 1088
#define IMAGE_HEIGHT 1288
#define DEGREE_TO_RADIAN(deg) ((deg) * M_PI / 180.0)

// 图像数据结构
typedef struct {
    int width;
    int height;
    unsigned char *data; // 指向图像像素数据的指针
} Image;

// 从RAW8文件读取图像数据
Image *read_raw8_image(const char *filename, int width, int height) {
    FILE *file = fopen(filename, "rb");
    if (!file) {
        printf("无法打开文件 %s\n", filename);
        return NULL;
    }

    // 在堆上为图像数据分配内存
    unsigned char *data = (unsigned char *)malloc(width * height * sizeof(unsigned char));
    if (!data) {
        printf("内存分配失败\n");
        fclose(file);
        return NULL;
    }

    // 读取图像数据
    fread(data, sizeof(unsigned char), width * height, file);
    fclose(file);

    // 创建Image结构并返回
    Image *image = (Image *)malloc(sizeof(Image));
    image->width = width;
    image->height = height;
    image->data = data;

    return image;
}

// 将图像数据写入RAW8文件
void write_raw8_image(const char *filename, Image *image) {
    FILE *file = fopen(filename, "wb");
    if (!file) {
        printf("无法打开文件 %s\n", filename);
        return;
    }

    // 写入图像数据
    fwrite(image->data, sizeof(unsigned char), image->width * image->height, file);

    fclose(file);
}
// 释放图像内存
void free_image(Image *image) {
    if (image) {
        if (image->data) {
            free(image->data);
        }
        free(image);
    }
}

// 仿射变换函数：旋转图像
Image *affine_transform(Image *image, float angle) {
    int width = image->width;
    int height = image->height;

    // 角度转弧度
    float radians = DEGREE_TO_RADIAN(angle);
    float cos_theta = cos(radians);
    float sin_theta = sin(radians);

    // 中心点
    int center_x = width / 2;
    int center_y = height / 2;

    // 在堆上为旋转后的图像数据分配内存
    Image *rotated_image = (Image *)malloc(sizeof(Image));
    rotated_image->width = width;
    rotated_image->height = height;
    rotated_image->data = (unsigned char *)malloc(width * height * sizeof(unsigned char));

    // 初始化旋转后图像为0（黑色）
    for (int i = 0; i < width * height; i++) {
        rotated_image->data[i] = 0;
    }

    // 遍历原始图像的每一个像素
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            // 计算相对于中心点的偏移
            int x_offset = x - center_x;
            int y_offset = y - center_y;

            // 进行仿射变换（旋转）
            int new_x = (int)(cos_theta * x_offset + sin_theta * y_offset) + center_x;
            int new_y = (int)(-sin_theta * x_offset + cos_theta * y_offset) + center_y;

            // 检查新坐标是否在图像范围内
            if (new_x >= 0 && new_x < width && new_y >= 0 && new_y < height) {
                rotated_image->data[new_y * width + new_x] = image->data[y * width + x];
            }
        }
    }

    return rotated_image;
}
int main() {
    // 读取RAW8图像
    Image *image = read_raw8_image("input.raw", IMAGE_WIDTH, IMAGE_HEIGHT);
    if (!image) {
        return 1;
    }

    // 对图像进行15度旋转
    Image *rotated_image = affine_transform(image, 15.0);

    // 保存旋转后的图像
    write_raw8_image("rotated_output.raw", rotated_image);

    // 释放内存
    free_image(image);
    free_image(rotated_image);

    return 0;
}

高性能处理之OpenCL实现

#include <stdio.h>
#include <stdlib.h>
#include <CL/cl.h>
#include <math.h>
#include <time.h>
#include <stdint.h>

#define CHECK_ERROR(err) if (err != CL_SUCCESS) { printf("OpenCL Error at line %d: %d\n", __LINE__, err); exit(EXIT_FAILURE); }

// OpenCL kernel source for affine transformation (rotation)
const char *kernel_source =
"__kernel void affine_rotate(__global const uchar *input, __global uchar *output, int width, int height, float cos_theta, float sin_theta, int new_width, int new_height) { \n"
"    int x = get_global_id(0); \n"
"    int y = get_global_id(1); \n"
"    float cx = width / 2.0f; \n"
"    float cy = height / 2.0f; \n"
"    float new_x = cos_theta * (x - new_width / 2.0f) + sin_theta * (y - new_height / 2.0f) + cx; \n"
"    float new_y = -sin_theta * (x - new_width / 2.0f) + cos_theta * (y - new_height / 2.0f) + cy; \n"
"    int src_x = (int)(new_x + 0.5f); \n"
"    int src_y = (int)(new_y + 0.5f); \n"
"    if (src_x >= 0 && src_x < width && src_y >= 0 && src_y < height) { \n"
"        output[y * new_width + x] = input[src_y * width + src_x]; \n"
"    } else { \n"
"        output[y * new_width + x] = 0; \n"
"    } \n"
"} \n";

static uint64_t get_monotonic_time_ns(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return (uint64_t)ts.tv_sec * 1000000000 + ts.tv_nsec;
}

void read_raw_image(const char *filename, unsigned char *buffer, size_t size) {
    FILE *file = fopen(filename, "rb");
    if (!file) {
        perror("Failed to open input file");
        exit(EXIT_FAILURE);
    }
    fread(buffer, 1, size, file);
    fclose(file);
}

void write_raw_image(const char *filename, const unsigned char *buffer, size_t size) {
    FILE *file = fopen(filename, "wb");
    if (!file) {
        perror("Failed to open output file");
        exit(EXIT_FAILURE);
    }
    fwrite(buffer, 1, size, file);
    fclose(file);
}

int main(int argc, char **argv) {
    uint64_t ts = 0;
    if (argc != 4) {
        printf("Usage: %s <input.raw> <output.raw> <rotation_angle_in_degrees>\n", argv[0]);
        return EXIT_FAILURE;
    }

    const char *input_filename = argv[1];
    const char *output_filename = argv[2];
    float angle_degrees = atof(argv[3]);
    float angle_radians = angle_degrees * M_PI / 180.0f;

    // 输入图像的宽度和高度
    const int width = 1088;
    const int height = 1288;
    const int new_width = 1088;  // 简单处理，这里不改变图片的尺寸
    const int new_height = 1288;

    // 分配输入和输出图像缓冲区
    size_t image_size = width * height;
    unsigned char *input_image = (unsigned char*)malloc(image_size);
    unsigned char *output_image = (unsigned char*)malloc(image_size);

    if (!input_image || !output_image) {
        printf("Failed to allocate memory for images.\n");
        return EXIT_FAILURE;
    }

    // 读取RAW8图像
    read_raw_image(input_filename, input_image, image_size);

    cl_platform_id platform;
    cl_device_id device;
    cl_context context;
    cl_command_queue queue;
    cl_program program;
    cl_kernel kernel;
    cl_mem buffer_input, buffer_output;
    cl_int err;

    // 初始化OpenCL
    err = clGetPlatformIDs(1, &platform, NULL);
    CHECK_ERROR(err);
    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
    CHECK_ERROR(err);
    context = clCreateContext(NULL, 1, &device, NULL, NULL, &err);
    CHECK_ERROR(err);
    queue = clCreateCommandQueue(context, device, CL_QUEUE_PROFILING_ENABLE, &err);
    CHECK_ERROR(err);

    // 创建内存缓冲区
    buffer_input = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, image_size, input_image, &err);
    CHECK_ERROR(err);
    buffer_output = clCreateBuffer(context, CL_MEM_WRITE_ONLY, image_size, NULL, &err);
    CHECK_ERROR(err);

    // 编译内核
    program = clCreateProgramWithSource(context, 1, &kernel_source, NULL, &err);
    CHECK_ERROR(err);
    err = clBuildProgram(program, 1, &device, NULL, NULL, NULL);
    if (err != CL_SUCCESS) {
        size_t log_size;
        clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
        char *log = (char *)malloc(log_size);
        clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, log, NULL);
        printf("Build log:\n%s\n", log);
        free(log);
        CHECK_ERROR(err);
    }

    // 创建内核
    kernel = clCreateKernel(program, "affine_rotate", &err);
    CHECK_ERROR(err);

    // 设置内核参数
    float cos_theta = cos(angle_radians);
    float sin_theta = sin(angle_radians);
    err = clSetKernelArg(kernel, 0, sizeof(cl_mem), &buffer_input);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 1, sizeof(cl_mem), &buffer_output);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 2, sizeof(int), &width);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 3, sizeof(int), &height);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 4, sizeof(float), &cos_theta);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 5, sizeof(float), &sin_theta);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 6, sizeof(int), &new_width);
    CHECK_ERROR(err);
    err = clSetKernelArg(kernel, 7, sizeof(int), &new_height);
    CHECK_ERROR(err);

    // 设置执行维度
    size_t global_size[2] = {new_width, new_height};

    ts =  get_monotonic_time_ns();
    // 记录执行时间
    cl_event event;
    err = clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global_size, NULL, 0, NULL, &event);
    CHECK_ERROR(err);

    // 等待完成
    err = clFinish(queue);
    CHECK_ERROR(err);

    // 获取执行时间
    cl_ulong time_start, time_end;
    err = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(time_start), &time_start, NULL);
    CHECK_ERROR(err);
    err = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(time_end), &time_end, NULL);
    CHECK_ERROR(err);

    double time_ns = (double)(time_end - time_start);
    printf("Rotation took %.3f ns\n", time_ns);

    // 读取处理后的图像
    err = clEnqueueReadBuffer(queue, buffer_output, CL_TRUE, 0, image_size, output_image, 0, NULL, NULL);
    CHECK_ERROR(err);
	printf("cost_ns:%llu\n", get_monotonic_time_ns() - ts);

    // 写入处理后的图像
    write_raw_image(output_filename, output_image, image_size);

    // 释放资源
    clReleaseMemObject(buffer_input);
    clReleaseMemObject(buffer_output);
    clReleaseProgram(program);
    clReleaseKernel(kernel);
    clReleaseCommandQueue(queue);
    clReleaseContext(context);
    clReleaseEvent(event);

    free(input_image);
    free(output_image);

    return 0;
}

在nxp 8plus的平台上处理1088x1288的raw8数据旋转15度大概耗时30ms;

优化思路之查表法(待解决问题）

但是目前下列代码还存在旋转以及图像问题，耗时也比上述的C代码在相同的nt98528平台上要大，上述C实验耗时大概170ms，这个查表法需要耗时250ms左右
问题1：当行列调转时，耗时从250ms降到了160ms，但是图像还存在问题；
问题2：图像显示问题；

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define ANGLE_DIVISIONS 1024  // 角度分割为1024份
#define MAX_SIDE_LENGTH 2048  // 最大边长支持2048个值
#define WIDTH 1088            // 图像宽度
#define HEIGHT 1288           // 图像高度

void rotate_image_using_table(unsigned char* src_image, unsigned char* dst_image, 
                              double **sin_table, double **cos_table, int angle_index) {
    int center_x = WIDTH / 2;
    int center_y = HEIGHT / 2;

    for (int y = 0; y < HEIGHT; y++) {
        for (int x = 0; x < WIDTH; x++) {
            // 计算相对于中心的坐标
            int trans_x = x - center_x;
            int trans_y = y - center_y;

            // 通过查表获取旋转后的坐标
            int new_x = (int)(cos_table[angle_index][abs(trans_x)] - sin_table[angle_index][abs(trans_y)]) + center_x;
            int new_y = (int)(cos_table[angle_index][abs(trans_y)] - sin_table[angle_index][abs(trans_x)]) + center_y;

            // 检查新坐标是否在图像边界内
            if (new_x >= 0 && new_x < WIDTH && new_y >= 0 && new_y < HEIGHT) {
                dst_image[new_y * WIDTH + new_x] = src_image[y * WIDTH + x];
            }
        }
    }
}

int main() {
    // 动态分配 sin 和 cos 表
    double **sin_table = (double **)malloc(ANGLE_DIVISIONS* sizeof(double *));
    double **cos_table = (double **)malloc(ANGLE_DIVISIONS* sizeof(double *));
    
    for (int i = 0; i < ANGLE_DIVISIONS; i++) {
        sin_table[i] = (double *)malloc(MAX_SIDE_LENGTH* sizeof(double));
        cos_table[i] = (double *)malloc(MAX_SIDE_LENGTH* sizeof(double));
    }

    // 计算每个坐标差值对应的 sin 和 cos 值，并存入表中
    double angle_increment = 90.0 / ANGLE_DIVISIONS;
    for (int i = 0; i < ANGLE_DIVISIONS; i++) {
        for (int j = 0; j < MAX_SIDE_LENGTH; j++) {
            double angle = i * angle_increment;
            sin_table[i][j] = j * sin(angle);
            cos_table[i][j] = j * cos(angle);
        }
    }

    // 假设我们要旋转15度
    int angle_index = (int)(15.0 / 90.0 * ANGLE_DIVISIONS); // 15度在1024份中的索引

    // 分配图像数据的内存
    unsigned char *src_image = (unsigned char *)malloc(WIDTH * HEIGHT * sizeof(unsigned char));
    unsigned char *dst_image = (unsigned char *)calloc(WIDTH * HEIGHT, sizeof(unsigned char));

    // 打开 input.raw 文件并读取数据
    FILE *input_file = fopen("input.raw", "rb");
    if (input_file == NULL) {
        perror("无法打开 input.raw 文件");
        return 1;
    }
    fread(src_image, sizeof(unsigned char), WIDTH * HEIGHT, input_file);
    fclose(input_file);

    // 统计旋转操作耗时
    struct timespec start_time, end_time;
    clock_gettime(CLOCK_MONOTONIC, &start_time);

    // 使用查表法旋转图像
    rotate_image_using_table(src_image, dst_image, sin_table, cos_table, angle_index);

    clock_gettime(CLOCK_MONOTONIC, &end_time);

    // 计算耗时
    long seconds = end_time.tv_sec - start_time.tv_sec;
    long nanoseconds = end_time.tv_nsec - start_time.tv_nsec;
    long total_time_ns = seconds * 1e9 + nanoseconds;
    printf("旋转操作耗时: %ld ns\n", total_time_ns);

    // 打开 deal.raw 文件并写入处理后的数据
    FILE *output_file = fopen("deal.raw", "wb");
    if (output_file == NULL) {
        perror("无法打开 deal.raw 文件");
        return 1;
    }
    fwrite(dst_image, sizeof(unsigned char), WIDTH * HEIGHT, output_file);
    fclose(output_file);

    // 释放内存
    free(src_image);
    free(dst_image);
    for (int i = 0; i < ANGLE_DIVISIONS; i++) {
        free(sin_table[i]);
        free(cos_table[i]);
    }
    free(sin_table);
    free(cos_table);

    return 0;
}