基于感知哈希的图像比较：用“纯C语言”编写。轶事观点
一篇轶事报告，描述了在已有35年历史的滚动开发中，使用纯C语言编写最小感知哈希功能时遇到的问题，目标是在新手Win 10/11安装上运行。
引言
首先：我真的指的是纯C语言，实际上是K&R C语言。为什么？后面会讨论。
其次：我对图像比较的理解不够深入，无法尝试解释图像比较的工作原理。
第三：我已经让图像比较（或者看起来像它的东西）工作起来了。（哦，天哪，这里有个糟糕的双关语）
第四：我将重点介绍我在寻找近似重复图像的过程中发现的一些结果。这些发现包括：
- 在纯C语言中使用
OleLoadPicturePath和IPicture很容易，但性能不太好。 - 在纯C语言中使用为C++和COM设计的
gdi+是可行的，但需要更多的工作，而且似乎比OleLoadPicturePath和IPicture性能更好。 - 离散余弦变换本质上是困难的，但在本文的背景下使用起来非常非常容易。
- 对于这个特定的用例，默认的离散余弦变换函数可以大大加速，达到150倍以上。
- 压缩位图时需要考虑半色调校正。
- 使用本文描述的方法，旋转后的图像不是原始图像的近似重复；-(
- 在我涉及30,000多张图像的特定案例中，数据库帮助巨大。
- 而Burkhard Keller树对我帮助更大。
在此过程中，你可能会注意到，关于纯C语言、逆波兰表示法后置Forth类语言以及另一个经典之作：Basic的组合对独立程序员的好处，虽然不是抱怨，但绝对是响亮的嘟囔。
虽然我怀疑我有超过10,000小时的编程经验，但我不是编程大师，所以本文中的所有内容都必须限定为据我所知（可能知道的并不多）。
链接和免责声明
我对感知哈希的理解来自在线搜索。以下是我用于获得有限理解的一些网站（受我能力限制，而非信息来源限制）
我与以下网站、其作者以及任何中间人或法律实体没有任何关联。我提供这份我访问过的一些网站列表，仅仅是为了在你想要快速入门相关网站时提供帮助。在撰写本文时，这些链接是有效的，我不知道浏览指定链接会产生任何恶意或有害影响。
- 432-Looks-Like-It
- how_to_automatically_identify_similar_images_using_phash
- phashtext
- open source perceptual hash library
- phash-determine-whether-2-images-are-equal-or-not
感谢与致谢
使用纯C语言是可以的，但在Windows环境下，你迟早会遇到使用组件对象模型（COM）。我使用纯C语言调用COM的入门知识来自CodeProject的文章。可以尝试Jeff Glatt于2006年撰写的COM-in-plain-C（不胜感激！）。
我使用gdi+的入门知识来自stackoverflow上的一个讨论，具体是：gdi+相关文章
我想借此机会赞美维护SQLite的人们。简单来说，这是一项了不起的工作。没有SQLite，我想我自2014年以来的职业生涯会少很多乐趣。我非常羡慕，也非常感激。
本文中我提到了一个元数据标签和浏览工具，我相信有很多这样的工具，但我有限的经验是使用exiftool。正是在创建我自己的元数据管理工具时，我偶然发现了exiftool。而正是这个元数据管理工具让我接触到了感知哈希。
为什么选择纯C语言
我对K&R C的增强功能没有任何问题，我想要表达的是：我不使用C++，无论是面向对象编程，还是（这是重要的部分）命令式编程。
为什么？因为C语言允许将指针值从一种指针类型复制到另一种指针类型，而无需任何显式类型转换表达式。由于我使用（稍后解释）一种Forth类虚拟机，我（可能很愚蠢地）大量使用了这种指针复制。这样的代码在C++环境中无法编译。我的错。但是，正如稍后在Direct2D问题中讨论的，我可能已经找到了一个平稳迁移到C++的途径。
我编程已经超过52年了，其中42年是专业编程。当为自己编程时，我在过程中获得巨大的乐趣，而在结果中获得的乐趣较少。用编程术语来说，这意味着理解和实现一个功能对我来说比使用最终程序更重要。如果我想要一个带有GUI的程序，我可以使用一个GUI框架，或者我可以在（Windows中）消息处理级别编写自己的GUI。猜猜我使用的是哪种方法！
我从1979年开始从事编程工作，使用接近硬件的语言：各种CPU的汇编语言（Intel i80/i85/i86..i486，Motorola 68000，Sun Sparc），以及“更高级”的接近硬件的语言，如C、PL/M和CHILL。
在此期间，我接触过（但未深入）：Fortran、Algol、Pascal、Smalltalk、C++（用于面向对象）和C#。1972年，Basic引起了我的注意，认为它是一种有趣的语言；1988年，Forth因其简洁性和潜力而令我震惊。老实说，我从未给Lisp机会，这是我的错，不是Lisp的错。对于Web编程，我会乐意使用PHP。Python、Haskell、Rust等等对我来说出现得太晚了。
我喜欢编写功能块，并拥有一个交互式环境，在那里我可以将这些块组合起来以解决问题。这很像最初的乐高概念。创建一组有限的基本构建块，只有当无法通过基本块的集合模拟时才创建自定义块。
自1988年以来，受Forth的影响，我一直遵循这个原则。我用C语言编写我的构建块，并将这些块放入一个单独的可执行文件中，该可执行文件包含一个类似于传统Forth系统的逆波兰表示法后缀解释器。然后，构建块可以通过Forth类包装器暴露给命令行解释器/用户。
我将整个概念/事物/想法/随你怎么说称为快速交互式编程环境技术，简称RipeTech，而当构建块和支持的Forth类代码编译后，就是我的rtForth可执行文件。
回归基础
大约在2014年，我参与了一个大型项目（数亿欧元）。这个项目不是公司日常业务的一部分。他们的IT部门完全忙于超过2000万活跃客户的日常业务。但是这个项目需要大量的临时数据分析。多年来我了解到，非日常业务项目提出的临时数据分析请求对于大多数IT部门来说是禁忌。
Excel文件（有些每周90万行800列）在项目周围飞来飞去，一些倒霉的个人正在编写VBA宏，试图从数据中获取意义。呃，这很痛苦。
由于公司政策/成本（以及成员的技能水平），我们无法安装Access。因此，为了改善我的工作条件，我使用了我rtForth中包含的SQLite源代码（单文件，纯C），并对所需的SQLite功能进行了一些包装。SQL意味着我可以大大简化和加快我的数据分析工作量。而SQLite意味着我可以进行无服务器操作。
但项目内的交换单位是Excel文件。所以我添加了C级别的组件对象模型（COM）功能，编写了一些rtForth包装器，并使用对象链接与嵌入/自动化/(COM)与Excel交换数据。
这对我有所帮助，但我相信99.9%的人在逆波兰表示法后缀编程方面存在问题。
于是我编写了一个脚本，当然是用rtForth，并包含一些用于提高性能的C语言构建块，这个脚本用一个基于Microsoft Visual Basic等Basic类解释器语法的内联解释器覆盖了rtForth。我将这些脚本嵌入到rtForth可执行文件中，最终得到了一个带有内置SQL和（COM）功能的独立Basic类解释器。这是一种具备Microsoft VBA技能的人可以轻松适应的语言。而且它是一个不需要安装的可执行文件，因此在公司PC/笔记本电脑上运行没有任何问题（事后来看，你可能会觉得这，呃，有点调皮！）。
我一个人，每年向项目交付多达1500次数据分析，这些分析基于每周更新的项目数据（90万行800列），并使用了SQLite/COM和脚本。我参与这个项目将近7年，非常愉快，依靠我的爱好过活。我真的生活在天堂里。
我称这个Basic类独立可执行文件为pwBasic，pw代表我用Basic编程的PlodWare。我早在1971/72年就写了我的第一个Basic程序，虽然我喜欢Basic，但我在摸索了40多年后才最终写出了自己的程序！
在我看来（我并非特别谦虚），这种功能/语言分层方法具有高度的灵活性。可以在适当的时候在C语言级别添加功能，然后用rtForth或pwBasic编写其包装器。
举例：我想为某个应用程序提供内联帮助，并希望相同的信息也能在线获取。因此我用HTML编写了帮助文本（一次），并将一个HTML浏览器控件作为一个构建块放入rtForth环境中。这个HTML块需要相当多的rtForth包装器，因为我希望提供对各种浏览器功能的交互式访问，例如设置URL、提取网页源代码、前进和后退浏览……现在，HTML构建块可用于未来的rtForth应用程序，通过添加pwBasic包装器，它也可以用于pwBasic应用程序。
当前的rtForth拥有超过1500个构建块，全部在一个可执行文件中，并且可编写脚本。
不，它不是一个巨大的可执行文件。小于3MB，单个文件，并且可以在新手Win 10或Win 11安装上运行。
即使是包含rtForth的pwBasic可执行文件也小于4MB。
结论
我用C语言编程，因为它迫使我审视细节问题。我避免使用框架，我承认如果目标是最终程序，框架非常棒，但如果目标是旅程，则不那么令人满意。
好吧，我不会尝试用自己的系统取代Windows操作系统，但我确实为多处理器平台编写过实时执行器。
而且我不会编写自己的代码将JPG文件渲染成位图，但这确实引起了我的兴趣，所以也许有一天我会踏上那段旅程。
我也用纯C语言编程，因为我在我的rtForth中投入了超过93,000行C语言源代码，这是许多许多快乐旅程的总和。
纯C语言限制了我吗？
有时会，但通常不会持续太久。
我曾想动态创建GUI，所以我添加了低级消息处理并编写了自己的控件处理（没有MFC或其他框架），我承认这并不容易，但却非常有用。
我想访问COM对象和COM服务器（例如Microsoft Excel），所以我编写了一个（出人意料地小巧的）构建块来使用CoCreateInstance和IDispatch等API和接口。
我甚至编写了一个迷你模块来调用.NET DLL/可执行文件中的代码。
正如我们将在本文后面看到的那样，我找到了一种在本质上是纯C语言的应用程序中混合使用C/C++的方法。
所以请不要因为我使用纯C语言或我的Forth和Basic版本而责备我。
不道歉。有时我觉得我使用这些语言是出于历史原因，因为我们这些老古董应该抱团取暖。
手套很合手！
图像比较：什么激发了我的兴趣
我曾经是一位热衷（尽管水平不高）的胶片摄影师，所以我有很多胶片图像（冲印照片、幻灯片、胶卷条：彩色和黑白，35毫米和6x6），当然还有一些从我的家人和我妻子的家人那里继承来的图像。总共有大约5000个胶片来源，再加上数码相机的图像，意味着我拥有超过30,000张图像。
作为一项退休活动，我决定将胶片图像数字化。这不是一项简单的任务，可能值得单独写一篇关于技巧和窍门的文章。
我认为为图像提供元数据会很好，这样我的孩子们就可以看着一张大约1909年的照片，并认出他们的曾曾祖母。
在纸质照片背面写字不是解决方案。最终，元数据可能会嵌入到数字文件中，作为EXIF或XMP标签（或其他）。
我决定，对于我所预期的需求，最灵活的解决方案是将元数据放入一个中央数据库，这样，当出现一个方便的图像文件标签解决方案时，我就可以简单地将元数据从数据库导出到标签工具，让它将标签插入到图像文件中。这个使用数据库的决定可能需要一些解释。简单来说，批量标签工具需要以特定格式（CSV、可能Excel、JSON或其他）提供标签信息。我尚不知道我需要什么格式，也不知道我可能需要哪些标签字段。因此，一个能够将其数据导出为CSV、Excel或JSON，或者为其创建自定义导出例程的数据库，对我来说似乎是一个明智的方法。
是的，我知道有元数据标记工具（例如exiftool等等），但目前，我希望数据库具有灵活性。那么，我所知道的哪个程序具有灵活的脚本功能和内置数据库呢？瞧，是的，pwBasic！
10,000行脚本之后，我拥有了一个灵活的伴随式元数据管理系统，并额外支持图像文件。这是元数据管理系统GUI界面之一的截图。
我用一个脚本（用pwBasic编写）完成了概念验证，将元数据从SQLite导出到CSV文件，然后触发exiftool将元数据写入图像文件中的标签字段。
冷酷的现实
遗憾的是，多年来，一些模拟图像被多次扫描/数字化，导致磁盘上存在几乎相同的副本。由于从旧（小）备份磁盘迁移到更大的集中存储系统，我还得到了某些图像文件的重复副本。
使用MD5校验和识别重复文件很容易。然而，几乎相同的图像（过去是，现在也是）有点棘手。
但大约三周前，我听说了一些数字图像比较技术，例如：
- 平均哈希
- 感知哈希
- 差异哈希
- 小波哈希
我对傅里叶变换、快速傅里叶变换和离散傅里叶变换的理解有限，所以我决定尝试使用离散余弦变换（DCT）的感知哈希技术。
不要被数学吓倒。DCT的编码确实很容易（源代码在下面提供），并且不需要任何参数调整。
但是在过去的三周里，我遇到了一些问题，想在这里分享一下。
创建感知哈希的步骤
首先，将模拟图像数字化。我在此不再赘述，它值得单独写一篇文章。
我将所有图像文件放在硬盘上的一个文件夹树中（为了性能，是一个固态设备），然后扫描该树以获取文件名，并将其放入SQLite数据库中。文件夹树的命名和结构可以简化元数据创建，这有一些好处，但这与本文无关。
然后我将数据库中引用（作为路径/名称）的每个图像文件提交到我pwBasic可执行文件中的一个C语言函数。该C函数执行以下操作：
- 将文件读取/渲染成位图（我尝试了
OleLoadPicturePath/IPICTURE API和gdi+ API） - 将位图压缩为32x32位图（使用
StretchBlt） - 将32x32位图从彩色转换为8位灰度（循环中的简单循环，但由于人类视觉感知而有趣）
- 对32x32灰度位图执行离散余弦变换，创建一个32x32的浮点值数组
- 从32x32 DCT值数组中提取一个特定的8x8子数组
- 获取8x8子数组的平均值
- 如果数组元素的值大于平均值，则将8x8数组中的值替换为
1，否则替换为0 - 创建一个64位值，其中第n位的值与8x8数组的第n个元素的值
1或0相同 - C函数然后返回64位值，即感知哈希
然后将64位值添加到相关图像的数据库条目中。
感知哈希创建完成。任务完成，呃，除了细节之外。
将图像文件读取/渲染到位图
感知哈希算法需要将图像视为灰度像素值的矩阵。但典型的图像文件（如JPG）将像素值保存在压缩、复杂的T数据结构中。幸运的是，在新手Windows系统中，有一些内置API可以将图像文件转换为简单的位图。当然，也有第三方库可以做到这一点，甚至更好！但我的目标是新手Windows系统。
最初，我使用OleLoadPicturePath读取图像文件以将其加载到IPicture结构中。然后使用IPicture的渲染方法将IPicture结构转换为位图。OleLoadPicturePath和IPicture是为COM对象设计的。但稍加努力，您就可以从纯C语言中调用它们。
// plain old C
// incomplete code, provided for guidelines only
// ------------------------------------------------------------------
// read the source IMAGE to a full size bitmap
// ------------------------------------------------------------------
// handle the filename
TCHAR * OlePathName ;
// Sorry, these two lines are rtForth code, they simply get the filename from a
// the underlying virtual machine.
fth_to_zero () ;
OlePathName = pop_ptr () ;
IPicture * Ipic = NULL;
HRESULT hr;
hr = OleLoadPicturePath
( OlePathName
, NULL
, 0
, 0
, & IID_IPicture
, (void *)(&Ipic)
) ;
if (hr)
{ // error handling, returns a ZERO to the underlying virtual machine
push_ucell (0) ;
return ;
}
// ------------------------------------------------------------------
HDC hDC = GetDC(NULL);
const int HIMETRIC_PER_INCH = 2540;
int nPixelsPerInchX = GetDeviceCaps(hDC, LOGPIXELSX);
int nPixelsPerInchY = GetDeviceCaps(hDC, LOGPIXELSY);
// get width and height of picture
SIZE sizeInHiMetric;
Ipic->lpVtbl->get_Width (Ipic,&sizeInHiMetric.cx);
Ipic->lpVtbl->get_Height(Ipic,&sizeInHiMetric.cy);
// convert himetric to pixels
myUInt targetWidth = (nPixelsPerInchX * sizeInHiMetric.cx +
HIMETRIC_PER_INCH / 2) / HIMETRIC_PER_INCH;
myUInt targetHeight = (nPixelsPerInchY * sizeInHiMetric.cy +
HIMETRIC_PER_INCH / 2) / HIMETRIC_PER_INCH;
// ------------------------------------------------------------------
// get a BITMAP big enough to hold the entire image
// Sorry I can not provide all of the code behind int_canvas_create_dib
canvasHANDLE * pDstCanvasHandle = int_canvas_create_dib
( targetWidth , targetHeight ) ;
if ( ! pDstCanvasHandle )
{ // error handling, release the IPicture ,
// then returns a ZERO to the underlying virtual machine
hr = Ipic->lpVtbl-> Release( Ipic ) ;
fth_zero () ;
return ;
}
// ------------------------------------------------------------------
// render image
hr = Ipic->lpVtbl->Render
( Ipic // the provider of the image pixel
, pDstCanvasHandle->dc // device context of the Bitmap to receive
// the rendered image
, 0 , 0
, targetWidth , targetHeight
, 0 // src x
, sizeInHiMetric.cy // src y
, sizeInHiMetric.cx
, -sizeInHiMetric.cy
, NULL
) ;
hr = Ipic->lpVtbl-> Release( Ipic ) ; // can release the IPicture,
// we only need the
// bitmap in pDstCanvasHandle
作为我另一个rtForth/pwBasic项目的一部分，我使用gdi+ API实现了一个图像读取器/渲染器，我相信它使用了硬件加速。
同样，这组API存在于COM环境中。虽然不是很困难，但它比OleLoadPicturePath和IPicture方法需要更多的代码。它需要：
GdiplusStartup：设置gdi+环境GdipCreateBitmapFromFile：将文件读取到gdi+位图GdipGetPropertyItemSize：获取gdi+位图的尺寸。用于传输到标准位图GdipCreateHBITMAPFromBitmap：将gdi+位图转换为标准位图GdipDisposeImage：释放gdi+图像gdipStop：释放gdi+环境
此外，我还遇到了一个问题，似乎在纯C语言中没有对gdi+的头文件支持。因此，我不得不查找/编写我需要的gdi+接口的函数原型，并手动编写一些数据结构和枚举。我从stackoverflow上的一个讨论中了解了如何使用gdi+，具体是：gdi+相关文章。
再次声明，不作任何道歉，这次是针对编码风格的。
这是我使用的gdi+代码。
// plain old C
#pragma comment(lib, "gdiplus.lib")
struct GdiplusStartupInputTYPE
{
UINT32 GdiplusVersion;
UINT_PTR DebugEventCallback;
BOOL SuppressBackgroundThread;
BOOL SuppressExternalCodecs;
VOID * fp ;
} ;
typedef struct GdiplusStartupInputTYPE GdiplusStartupInput ;
typedef struct
{
DWORD id;
ULONG length;
WORD type;
VOID* value;
} PropertyItem;
typedef enum
{
RotateNoneFlipNone = 0,
Rotate90FlipNone = 1,
Rotate180FlipNone = 2,
Rotate270FlipNone = 3,
RotateNoneFlipX = 4,
Rotate90FlipX = 5,
Rotate180FlipX = 6,
Rotate270FlipX = 7,
RotateNoneFlipY = Rotate180FlipX,
Rotate90FlipY = Rotate270FlipX,
Rotate180FlipY = RotateNoneFlipX,
Rotate270FlipY = Rotate90FlipX,
RotateNoneFlipXY = Rotate180FlipNone,
Rotate90FlipXY = Rotate270FlipNone,
Rotate180FlipXY = RotateNoneFlipNone,
Rotate270FlipXY = Rotate90FlipNone
} RotateFlipType;
// in x64 "C" and CPP have the same
// calling convention so __stdcall is not needed (as it is in x32)
int GdiplusStartup(ULONG_PTR *token,
struct GdiplusStartupInput *gstart, struct GdiplusStartupInput *gstartout);
int GdiplusShutdown(ULONG_PTR token);
int GdipSaveImageToFile
(void *GpBitmap, WCHAR *filename, CLSID *encoder, void *params);
int GdipDisposeImage(void *GpBitmap);
int GdipCreateBitmapFromHBITMAP(HBITMAP bm, HPALETTE hpal, void *GpBitmap);
int GdipCreateBitmapFromFile(WCHAR *filename, void **GpBitmap);
int GdipGetPropertyItem (void *GpBitmap, PROPID propID, UINT size, void *buffer);
int GdipGetPropertyItemSize (void *GpBitmap, PROPID propID, UINT *size);
int GdipImageRotateFlip (void *GpBitmap , RotateFlipType type) ;
ULONG_PTR gdiplusToken;
uint64_t gdipStatus = 0 ;
uint64_t gdipStart ()
{
if ( gdipStatus ) return gdipStatus;
// get GDI+ going
GdiplusStartupInput gdiplusStartupInput;
gdiplusStartupInput.GdiplusVersion = 1;
gdiplusStartupInput.DebugEventCallback = NULL;
gdiplusStartupInput.SuppressBackgroundThread = 0;
gdiplusStartupInput.SuppressExternalCodecs = 0;
gdiplusStartupInput. fp = 0;
if ( GdiplusStartup(&gdiplusToken, &gdiplusStartupInput, NULL) )
{
return 0 ;
}
gdipStatus = 1;
return gdipStatus ;
}
void gdipStop ()
{
static uint64_t status = 0 ;
if ( ! gdipStatus ) return ;
GdiplusShutdown(gdiplusToken);
gdipStatus = 0;
}
void int_create_bitmap_from_file_GDIP ()
{
/*
** this "C" function is embedded in the "Forth" like virtual machine
** It obtains its input arguments by popping them off the vitual machine stack
** It returns values by pushing them on to the virtual machine stack
** The code is written in 64 bit mode, with wide characters as default
*/
// ------------------------------------------------------------------
// handle the filename
TCHAR * FileName ;
// sorry, these two are calls to the "Forth" like virtual machine
// they simply get the user provided image filename
// as a valid "C" string (zero terminated)
fth_to_zero () ;
FileName = pop_ptr () ;
int retStatus;
void * GpBitmap = NULL ;
// ------------------------------------------------------------------
// get GDI+ going
if ( ! gdipStart () )
{
fth_zero () ;
return ;
}
// ------------------------------------------------------------------
// get FILE
retStatus = GdipCreateBitmapFromFile ( FileName , & GpBitmap ) ;
// ------------------------------------------------------------------
// get FILE properties:
PropertyItem * itemProperties ;
myUInt wantedProperty ;
myUInt propertySize ;
// ------------------------------------------------------------------
// ORIENTATION
unsigned short ExifOrientation;
ExifOrientation = 1 ; // default
unsigned short gdipOrientation = 0 ; // default
wantedProperty = 0x0112 ; // the property id for the EXFI orientation meta-data
propertySize = 0 ;
retStatus = GdipGetPropertyItemSize ( GpBitmap, wantedProperty, & propertySize );
if ( retStatus EQU 0 )
{
itemProperties = (PropertyItem *)malloc(propertySize);
retStatus = GdipGetPropertyItem
( GpBitmap, wantedProperty, propertySize, itemProperties);
ExifOrientation = *( unsigned short * ) ( itemProperties->value ) ;
free ( itemProperties ) ;
switch (ExifOrientation )
{
case 1: // none
gdipOrientation = RotateNoneFlipNone ;
break ;
case 2: // flip horizontal
gdipOrientation = RotateNoneFlipX ;
break ;
case 3: // rotate 180
gdipOrientation = Rotate180FlipNone ;
break ;
case 4: // flip vertical
gdipOrientation = Rotate180FlipX ;
break ;
case 5: // transpos
gdipOrientation = Rotate90FlipX ;
break ;
case 6: // rotate 90
gdipOrientation = Rotate90FlipNone ;
break ;
case 7: // transverse
gdipOrientation = Rotate270FlipX ;
break ;
case 8: // rotate 270
gdipOrientation = Rotate270FlipNone ;
break ;
}
retStatus = GdipImageRotateFlip (GpBitmap , gdipOrientation) ;
}
// ------------------------------------------------------------------
// get FILE as HBITMAP
HBITMAP hBitmap = NULL;
retStatus = 0 ;
if ( retStatus EQU 0 )
{
retStatus = GdipCreateHBITMAPFromBitmap( GpBitmap, & hBitmap, 0) ;
}
if ( GpBitmap NEQU NULL )
GdipDisposeImage(GpBitmap); // destroy the 'Image' object
gdipStop () ;
// ------------------------------------------------------------------
// the only variable that we carry forward is hBitmap
// ------------------------------------------------------------------
push_ptr (hBitmap ); // pass to caller via virtual machine stack
}
令人欣慰的是，gdi+应该使用一些硬件加速。我不知道它是否使用了，但在我的笔记本电脑上测试我的整个图像文件集合时，gdi+似乎比OleLoadPicturePath和IPicture方法快了将近一倍。
图像压缩，或位图与半色调的问题
抱歉，如前所述，我无法解释感知哈希的工作原理，但我读过可以将图像压缩成32x32灰度像素图像，然后对32x32灰度图像进行DCT。因为DCT本质上是矩阵计算，所以对32x32矩阵进行DCT比对原始图像的位图矩阵（很容易是数千x数千）进行DCT的计算成本更低。但是，当我们减小矩阵大小时，我们会失去精度/清晰度。“他们”，无论“他们”是谁，都说32x32是可以的！
正是在这一点上，我确实遇到了问题。
当对照我认为近似相同的图像检查我的实现时，以下图像对的相似度很低。我当时想，它们肯定会有非常相似的感知哈希值。错了！我只测量到64%的相似度，也就是说，它们的感知哈希值在64位中相差23位。
这让我认为我的代码出了严重问题。我撞了无数次“南墙”。我睡得很不好。然后我编写了一个pwBasic脚本来并排比较一些位图。接着我将步骤分解：读取/渲染、压缩、灰度、DCT成块，并在过程中的各个点显示位图。
如果你将图像文件读取/渲染成位图（例如使用前面介绍的int_create_bitmap_from_file_GDIP函数），然后简单地将位图扔到屏幕上，如果位图需要缩小以适应屏幕，它将显得有些奇怪。
用Windows技术术语来说：通过StretchBlt将源位图绘制到WM_DRAWITEM消息提供的设备上下文中来响应WM_DRAWITEM消息。
// plain old C
rCode = StretchBlt
( destinationDeviceContext
, 0
, 0
, 32
, 32
, sourceDeviceContext
, 0
, 0
, widthOfSourceBitmap
, heigthOfSourceBitmap
, SRCCOPY
) ;
对于城堡的左图（上图）：左边（下图）是你（可能）想要的，但右边是你简单地把位图扔到屏幕上的结果！
对于城堡的右图：左边（下图）是您（可能）想要的，但右边是您直接将位图扔到屏幕上的结果！
那又怎样，谁在乎
好吧，结果我关心，我想！
我最终意识到，从全尺寸图像位图到DCT计算所需的32x32位图的压缩，没有处理半色调。
（下图）是（未经半色调处理）代码链正在比较的内容，仅显示64%的相似度，即它们的感知哈希值在64位中相差23位。
当我使用半色调校正（如果这个说法正确的话？）时，我得到了更好的79%相似度，也就是说，它们的感知哈希值在64位中相差12位。
带有半色调校正的城堡图像比较如下图所示。
带或不带半色调校正
英国有句俗话说“一燕不成夏”（燕子：一种候鸟：“英国”的燕子在非洲过冬，夏天返回英国）。
所以我测试了许多带有和不带有半色调校正的图像对。总的来说，在我看来使用半色调校正更好。
所需要做的就是在调用StretchBlt之前立即添加一行代码SetStretchBltMode，该StretchBlt将全尺寸图像位图压缩到用于离散余弦变换的32x32像素位图。
实际上，Microsoft文档说要添加两行SetStretchBltMode和SetBrushOrgEx。
// plain old C
SetStretchBltMode(destinationDeviceContext HALFTONE); // correction for halftone
// issue, HALFTONE = 0x0004
SetBrushOrgEx(destinationDeviceContext , 0, 0, NULL); // recommended by Microsoft
// documentation!
rCode = StretchBlt
( destinationDeviceContext
, 0
, 0
, 32
, 32
, sourceDeviceContext
, 0
, 0
, widthOfSourceBitmap
, heigthOfSourceBitmap
, SRCCOPY
) ;
扭转局面
后来，在我开发这种方法的三周时间里，我注意到一些近似相同的图像具有截然不同的感知哈希值，这主要是因为图像以不同的方向扫描。扫描图像的“方向”标签（如果它们有这样的标签的话）值不正确，因为它们无法知道图像的方向。
不知何故，我曾认为DCT算法会与方向无关。显然，我错了！
这意味着图像方向与基于感知哈希的图像比较相关（至少在我的代码实现中如此）。
为了好玩/参考，我展示了第一张城堡图像与它本身顺时针旋转90度的版本进行比较，它们的感知哈希值在64位中相差36位，因此它们的相似度为36%。
所以，我决定添加一些代码来纠正图像方向。
gdi+方法需要另外两个gdi+方法（它们在上面的代码片段中）
GdipGetPropertyItem：获取图像方向GdipImageRotateFlip：将图像“扶正”
OleLoadPicturePath和IPicture的方法需要更多一点的代码（这是一种英式低调的说法！）。
- 如果您的编程环境支持，请对Shell.Application的
GetDetailsOfAPI进行脚本调用。 - 或者，如果你的编程环境支持，请调用Windows图像采集的脚本接口
- 或使用Windows图像组件查询典型JPG文件中嵌入的元数据属性，从而获取方向。
WIC是为COM接口设计的，在纯C语言中这是一个小问题。不管怎样，我编写了这个解决方案。
然后，一旦你获得了方向，就相应地旋转图像。旋转代码并不困难，但我花了相当长的时间才把标记为???的行中的逻辑搞对，见下文。
// plain old C
uTYPE32 * srcPtr , * srcDibData = pSrcCanvasHandle->data; // ptr to the pixel
// data in the source
// bitmap
uTYPE32 * dstPtr , * dstDibData = pCanvasHandleDst->data; // ptr to the pixel
// data in the
// destination bitmap
int srcStride = pSrcCanvasHandle->stride / sizeof(uTYPE32) ; // distance between
// two consecutive
// source pixel rows,
// in 32bit steps
int dstStride = pCanvasHandleDst->stride / sizeof(uTYPE32) ; // distance between
// two consecutive
// destination pixel
// rows, in 32bit steps
int row, col
srcPtr = srcDibData ;
dstPtr = dstDibData ;
switch ( rotation )
{
case 0:
for ( row = 0 ; row < pSrcCanvasHandle->max_h ; row++ )
{
for ( col = 0 ; col < pSrcCanvasHandle->max_w ; col++ )
{
*dstPtr++ = *srcPtr++ ; // ??? easy
}
}
break ;
case 90:
for ( row = 0 ; row < pSrcCanvasHandle->max_h ; row++ )
{
for ( col = 0 ; col < pSrcCanvasHandle->max_w ; col++ )
{
dstPtr [ (dstStride*(pSrcCanvasHandle->max_w - col-1) )
+ (row) ] = *srcPtr++ ; // ???
}
}
break;
case 180:
for ( row = 0 ; row < pSrcCanvasHandle->max_h ; row++ )
{
dstPtr = dstDibData +
(pSrcCanvasHandle->max_h - row -0)*dstStride ; // ???
for ( col = 0 ; col < pSrcCanvasHandle->max_w ; col++ )
{
*--dstPtr = *srcPtr++ ; // ???
}
}
break ;
case 270:
for ( row = 0 ; row < pSrcCanvasHandle->max_h ; row++ )
{
for ( col = 0 ; col < pSrcCanvasHandle->max_w ; col++ )
{
dstPtr [ (dstStride*col ) +
(pSrcCanvasHandle->max_h - row -1) ] = *srcPtr++ ; // ???
}
}
break;
default:
break;
}
或者，这就是巧妙之处。忘记使用OleLoadPicturePath和IPicture的方向校正方法，坚持使用gdi+方法，它只是更快（也更容易）。
将压缩的RGB位图转换为灰度
有一件事总是让我有点烦恼，那就是英式英语或美式英语拼写。我出生于英国，拥有双重国籍（德国）。但坦率地说，我的头发是grey，而不是gray。这对于一个年长的人来说是一种自然色（不是color）（我于1971年用CESIL编程语言写了我的第一个程序）。
Grey或gray，随便。我用于灰度转换的代码是
// plain old C
unsigned char * dibData = pReducedCanvasHandle->data; // this is a ptr to the
// pixel values in the
// 32 by 32 compressed
// bitmap
unsigned char * dptr ;
double grayArray [ 32 * 32 ] ;
double * dblPtr ;
dblPtr = grayArray ;
for ( row = 0 ; row < 32 ; row++ )
{
dptr = dibData + ( row * ( pReducedCanvasHandle->stride ) ) ;// this is a ptr
// to the 32 by 32
// compressed
// bitmap STRIDE
// value
// the distance
// in bytes between
// 2 consecutive
// rows of pixel
// data
for ( col = 0; col < 32 ; col ++ )
{
// the grayscale value is a weighted sum of the Red Green Blue components
// why these weights? Do not ask me, they are related to "human perception".
// And I did say it once or twice earlier on,
// do not ask me for explanations ;-)
*dblPtr++ = ( 0.072 * dptr[0] ) +
( 0.715 * dptr[1] ) +
( 0.213 * dptr[3] ) ;
dptr += sizeof(uTYPE32) ;
}
}
离散余弦变换
我将不尝试解释DCT的原理，原因很简单。我的数字信号处理知识在过去的42年中已经生疏到毫无价值的地步。十五年前，我曾尝试过离散傅里叶变换，只是为了好玩，创建了一个音频信号的瀑布显示，并修改了一个Goertzel滤波器作为摩尔斯电码解码器的一部分。但抱歉，我甚至无法开始解释其中的任何内容！
使用DCT比我预期的要容易得多。DCT的代码很容易在网上找到，所以我复制了一份并对其进行了一些优化。
// plain old C
double dct (double *DCTMatrix, double *ImgMatrix, int N, int M)
{
// This is basically what one finds online
// for a Discrete Cosine Transform algorithm in 2 dimensions.
// Oh, not strictly true, I have included code to find the
// average value of the DCTMatrix, excluding the value at DCTMatrix[0],
// apparently the value at DCTMatrix[0] is not helpful?
// The DCTMatrix and ImgMatrix have the same size N * M
// This code calls the innermost loop 1048576 times when N = M = 32
// On my laptop 1000 iterations (call to this function with image data)
// took 12,706 ms i.e. almost 13 seconds
int i, j, u, v;
int cnt = 0 ;
double DCTsum = 0.0 ;
for (u = 0; u < N; ++u)
{
for (v = 0; v < M; ++v)
{
DCTMatrix[(u*N) + v] = 0;
for (i = 0; i < N; i++)
{
for (j = 0; j < M; j++)
{
DCTMatrix[(u*N) + v] += ImgMatrix[(i*N) + j]
* cos(M_PI/((double)N)*(i+1./2.)*u)
* cos(M_PI/((double)M)*(j+1./2.)*v);
cnt=cnt++ ;
}
}
DCTsum += DCTMatrix[(u*N) + v] ;
}
}
DCTsum -= DCTMatrix[0] ;
return DCTsum / cnt ;
}
我对它进行了图像数据性能测试，1000次迭代耗时12.7秒。每张图像都需要执行DCT循环1048576次（32*32*32*32）。
还有很大的改进空间：无论是在可理解性还是性能方面。极其相关的是，虽然DCT是32x32，但我们只需要8x8子矩阵（行0:7，列0:7）。还可以通过简单的查找预定义表来替换昂贵的cos（余弦）调用。
经过几次优化尝试，我得到了以下结果：
// plain old C
int CoeffFlag = 0 ;
double Coeff [32][32] ;
double dctOptimised (double *DCTMatrix, double *IMGMatrix )
{
/*
Concerning: the expression
dctRCvalue+= IMGMatrix[(imgRow * RowMax) + imgCol] // this is dependant
// on the IMAGE
* cos( PIover32 *(imgRow+0.5)*dctRow) // this is a set of FIXED values
* cos( PIover32 *(imgCol+0.5)*dctCol); // this is a set of FIXED values
Let us call the 2 sets of FIXED values rCoeff and cCoeff
they both have the same set of values
= cos ( PIover32 * ( x + 0.5 ) * y ) for x = 0 .. 31 and y = for x = 0 .. 31
= 32*32 distinct COSINE values
= we could calculate these COSINE values in advance,
place them in an array, and (hopefully)
speed things up by doing a simple look up .
*/
#define PIover32 0.09817477042
#define RowMax 32
#define ColMax 32
int imgRow, imgCol ;
int dctRow, dctCol ;
int x , y ;
int cnt = 0 ;
double DCTsum = 0.0 ;
double dctRCvalue = 0.0 ;
if (! CoeffFlag )
{
for (x = 0 ; x < 32 ; x++ )
{
for (y = 0 ; y < 32 ;y++)
{
Coeff[x][y] = cos ( PIover32 * ( x + 0.5 ) * y ) ;
}
}
CoeffFlag = 1;
}
for (dctRow = 0; dctRow < 8; dctRow++)
{
for (dctCol = 0; dctCol < 8; dctCol++)
{
dctRCvalue = 0;
for (imgRow = 0; imgRow < RowMax; imgRow++)
{
for (imgCol = 0; imgCol < ColMax; imgCol++)
{
dctRCvalue+= IMGMatrix[(imgRow * RowMax) + imgCol]
* Coeff[imgRow][dctRow] // cos( PIover32 *(imgRow+0.5)*dctRow)
* Coeff[imgCol][dctCol] ; // cos( PIover32 *(imgCol+0.5)*dctCol) ;
cnt=cnt++ ;
}
}
DCTMatrix[(dctRow*RowMax) + dctCol] = dctRCvalue;
DCTsum += dctRCvalue;
}
}
DCTsum -= DCTMatrix[0] ;
return DCTsum / cnt ;
}
在我的笔记本电脑上，优化后的DCT进行1000次迭代耗时71毫秒，即0.071秒。与原始版本相比，这快了大约150倍。因此，30,000次未优化的DCT将耗时约390秒（5.5分钟），而优化后的DCT仅需21秒。每张图像都需要执行DCT循环65536次（32*32*8*8）。
所以为30,000多张图像节省了5分多钟，虽然不多，但每一点帮助都值得！
这是个好消息，但真正的开销似乎是读取/渲染图像文件到位图。我没有进行任何测量来验证这个假设。
Direct2D带来的读取/渲染性能提升
好吧，使用Direct2D可能会有性能提升。也许吧！我不知道（目前）。
我已经研究过使用Direct2D API，这些API应该很快。但是最新的Microsoft Windows SDK有一个不幸的特点，那就是Direct2D头文件不支持纯C语言。而且似乎Direct2D的一些接口方法很难映射到纯C语言。我搞不懂为什么会这样，但它就是这样！解决办法是在C++环境中编写调用Direct2D的代码。它不严格需要面向对象和类，只需要COM。
我害怕将我的代码从纯C语言移植到C++。至少因为我在93,000多行纯C语言源代码中大约有1500个函数调用需要在C++中编译之前进行处理。
在线搜索显示你可以混合使用纯C语言和C++，但是
- 跨C/C++边界调用函数需要在函数原型中格外小心。
WinMain代码应该用C++实现，以确保链接到正确的C++库。
嗯，不一定如此！至少如果你只需要C++和COM，而不需要对象类的话。
我用C++编写了我的Direct2D代码，并通过一个简单的构造让C++代码看到rtForth其余部分的C语言头文件。
// C++
// this is the top few lines of the C++ code
#include "f_d2d1.hpp" // C++ header file for the Direct2D code
#include < iostream> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < stdio.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < stdlib.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < string.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < Windows.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < d2d1.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
#include < wincodec.h> // SORRY, remove the space after the less-than sign.
// Seems to be a HTML issue
extern "C" {
#include "f_protos.h" // header file for those plain "C" functions that
// the Direct2D code might call
// no alteration to this header where needed
}
#pragma comment(lib, "d2d1.lib")
#pragma comment(lib, "Windowscodecs.lib")
将调用Direct2D的功能写入一个C++文件。然后，当一个纯C语言函数需要调用C++文件代码中的函数时，你需要在纯C语言源文件中包含一个头文件，其中包含所需C++函数的函数原型。例如：
// plain old C
#include "f_d2d1.hpp" // "C" header file for "C" source files,
// the header references the needed "C++" functions.
在f_d2d1.hpp文件中，“C++”函数原型通过以下方式提供：
// plain old C
#pragma once
#ifdef __cplusplus
extern "C" {
#endif
void fth_d2d1_main(); // functional prototype of a C++ function
// that is callable from plain old "C"
#ifdef __cplusplus
}
#endif
收尾工作
大胆地说，我可以在下面的代码片段之后结束。但是，还有一些后续问题。
// plain old C
int64_t hash_matrix ( double * matrixPtr, double averageValue )
{
int row, col ;
int64_t phash ;
for ( row = 0 ; row < 8 ; row++ )
{
for ( col = 0; col < 8 ; col++ )
{
if ( averageValue < matrixPtr[(row*32) + col ] )
phash = (phash<<1) + 1;
else
phash = phash<<1;
}
}
return phash ;
}
// and THEN, in some function that manages the steps,
// I have the following to finalise the perceptual hash generation
double avg8x8 ;
double * DCTMatrix [32 * 32 ] ; // matrix that computeDTC
// will fill with DCT calculated values
avg8x8 = computeDTC ( DCTMatrix, grayArray , 32, 32 ) ; // calculate the DCT
// from the grayscale
// bitmap
return hash_matrix ( DCTMatrix, avg8x8 ) ) ; // convert the DCT matrix to
// a 64 bit set of 1s and 0s,
// the perceptual hash
后续问题
定义“相似”
如果你想知道两张图片有多相似，那么获取它们的64位感知哈希值，对这两个值进行异或位运算，然后简单地计算异或结果中置位的位数（称为汉明距离）。我的经验（ admittedly limited to the last few weeks）是，汉明距离为12或更少表示相似！（如本文前面城堡比较所示。）
如果你搜索caltech 101 images dataset，你应该能找到一个包含大约9000张（小）图像的zip文件，这些图像具有不同程度的相似性。这里是它的链接：Caltech 101 images dataset。我下载了他们的图像，并编写了一个pwBasic脚本（不到200行），这样我就可以将图像拖放到屏幕上的两个区域之一，然后为这两张图像创建感知哈希。显示感知哈希为8x8矩阵，显示两张图像之间的汉明距离，并将该汉明距离转换为百分比相似度。
让我有点难过的是，我没有在网上找到一套完整的参考图像，包含它们的感知哈希值。我对自己工作的信心不足，希望能够对照公认的参考进行检查。可能我误解了Caltech 101 images dataset提供的内容，因为除了图像之外，zip文件还包含MATLAB文件中的信息，对我来说毫无意义。
考虑使用布克哈德-凯勒树
即使您的图像已在数据库中被引用，并带有感知哈希值，问题也并未结束。如果您要求一个现成的关系数据库管理系统（RDBMS）为其30,000多个条目中的每一个找到最近的邻居，其性能可能不会太好。
所以，我为我在pwBasic中使用的SQLite数据库实现了一个Burkhard Keller树扩展。尽管插入BK树值的代码不大（约60行纯C语言代码），并且查找最近邻居的代码也只有约80行，但要获得正确的逻辑却非常痛苦。我可能会就此写一篇文章（我想我甚至可以解释BK树的工作原理）。
但实现BK树的挫折是值得的。
性能：感知哈希和BK树
扫描30000多张图像，总计约82GB，以创建一个包含文件路径/名称和每个文件的感知哈希值的数据库，耗时61分钟（我的笔记本电脑是6年前的i7-7600U CPU @ 2.80GHz，代码未采用多线程）。然后，调用SQLite的Burkhard Keller树扩展以获取每张/所有图像的最近邻居，又耗时30秒。
纯粹作为观察。我儿子有一台全新的（1周大的）笔记本电脑，配备AMD Ryzen Pro 7 CPU，我们把图像文件和pwBasic可执行文件以及脚本复制（无需安装）到他的笔记本电脑上并重新运行了测试。必须强调的是，pwBasic是为单核（即单线程）编写的，所以我们没有期望太多的处理性能提升，两台笔记本电脑都配备了NVME固态硬盘，所以我们没有期望任何磁盘访问优势。我儿子的笔记本电脑耗时44分钟（我的耗时61分钟）。我儿子此后说服我认真考虑生成一些线程，将图像渲染负载分散到CPU核心上。（这对我来说听起来有点现代。）
后续步骤
我正在将感知哈希代码和Burkhard Keller树代码集成到我的pwBasic元数据管理系统脚本中，纯C语言代码已经在rtForth中并由pwBasic包装，我只需要决定如何在脚本中最好地使用它。
稍后，当pwBasic元数据管理系统真正稳定后，我将尝试Direct2D和C++，甚至可能将我纯C语言代码的其他部分移植到C++。虽然我发现当前的感知哈希计算性能可以接受，但为了计划中的元数据管理系统增强功能，我希望获得性能更好的图像渲染。
哦，是的，我会考虑多线程以及它是否/如何适应rtForth / pwBasic，但我对此存疑。
我将把SQLite Burkhard Keller树扩展源代码发布到公共领域，最初发布在我的网站上（https://www.ripetech.com）或（https://www.plodware.com，技术上是RipeTech的子域名）。
当我满意时，我将在我的网站上发布pwBasic元数据管理系统。
序言
正如作家道格拉斯·亚当斯在他精彩的著作《银河系漫游指南》中所写，天狼星控制论公司投诉部有以下非常恰当（？）的座右铭。
分享与享受
我想引用一句中国古语：与其诅咒黑暗，不如点亮蜡烛
我希望我的蜡烛，这篇文章，能带来一丝光明。
谨此谦卑。
GeorgeDS
历史
- 2023年6月4日：首次发布