PHASH: what changed

Revision 0 to 1. Removed lines have a minus sign; added lines have a plus sign.

Title: PHASH

New		Line
1	+	## Calculating Perceptual hashes (PHASH) and comparing image similarity
2	+
3	+	Here is a program to calculate the perceptual hash for images different than a CRC or md5 hash a perceptual hash will be similar for similar images. A slight change in the image (sharpen, blur, resize) will give perceptual hash that is similar where a CRC or MD5 hash will be completely different. I looked for a PHASH program for pascal and could not find any. I did find three methods for getting a perceptual hash. The first was using a DCT transform but the code was written in C and was not easy to link the program in. The second way is to resize it and make a hash based on whether each pixel was more or less then its next neighbor. The third way was to use image moments. ImageMagick(IM) uses this method but I found that after a certain threshold the images where radically different. Also calculating the differences required a lot of computing power as you had to square the difference of 42 numbers for each image comparison. I tried all three methods, only the first had good promise. The last one was not even part of the IM Lazarus interface.
4	+
5	+	This program uses ImageMagick and I want to express my appreciation for the work of the translators of this tool to work with Lazarus. However I had a lot of problems getting it to work properly. First because it has not kept up with the code there were some routines that did not work because they have been deprecated (and no longer available) in ImageMagick. Second because some of the newer routines were missing entirely. And third because I was using the Q8 version of ImageMagick and it appeared that the translation was for the Q16 version. PPixelPacket for instance assumes the colors are 16 bit and not 8 bit so I had to change how I handled the structure). I use 8 bit ImageMagick because 99.99999% of the pictures I wanted to process are 8 bit or less and the Q8 program is much faster then the Q16 version.
6	+
7	+	To get this to work under windows 7 you need to compile under 32bit mode (and call 32 bit ImageMagick) as the calling conventions for 64 bit are quite different and Lazarus does not support the 64 bit calling conventions. (this may change but keep it in mind).
8	+
9	+	This program takes one main argument - the directory to parse for images (currently only JPG and GIF but you can add more). The program maintains a file 'yalls' that contains information for each image it finds in that directory - deleting files that no longer exists, adding new files and if the md5 does not match recalculating data for that file. (for speed you can have it skip the md5 check and just rewrite existing lines for files by placing SKIP before the directory argument. Since I had a specific set of directories to parse if you do not use any argument it will parse those directories. The output file will consist of one line for each image it finds consisting of fields separated by blanks. The program expects the filenames are without spaces and would need modification otherwise. There was a memory leak with the program and after about 200,000 images it crashed. This may be gone but if the program dies just move the yall file to yalls and restart with JPEGINFO SKIP to restart from where it ended.
10	+
11	+	The output file contains:
12	+
13	+	` MD5 - 32 hex bytes of the MD5 sum of the entire file.`\
14	+	` Signature - 64 hex bytes that comes from IM identify that is I think the md5 of the image only`\
15	+	` PHASH - 16 hex bytes of the DCT hash of the image`\
16	+	` Date of the file in YYYYMMDDHHMMSS format (from ImageMagick identify)`\
17	+	` Size in bytes of the file`\
18	+	` Bitdepth of the image`\
19	+	` Width of the image`\
20	+	` Height of the image`\
21	+	` Filename of the image`\
22	+	` Complete fileid of the image`\
23	+	` Example:`\
24	+	` md5 fb54e12c44d2b14d4d3266251109799e`\
25	+	` signature 2318170444932902d0541895584df4c2f8f38ad8b1d4dad899870760817422f2`\
26	+	` phash CCFB4EAEF3F30C0D`\
27	+	` date 20151001010847`\
28	+	` size 901`\
29	+	` bitdepth 8`\
30	+	` width 32`\
31	+	` width 32`\
32	+	` filename 01.gif`\
33	+	` fid F:\CCC\01.gif`
34	+
35	+	I will add another program later if there is interest that calculates phash similarity. This program was written on a windows7 platform. Basically to get phash similarity you count the number of bits that are different. If that number is 0 the pictures are likely to be the same, if it is over 10 the pictures are likely to be different, otherwise the pictures may be similar to some degree. Note it takes about 24 hours to process 700000 images with this program on my computer. I estimate it may take another 24 hours to compare these hashes up to 4 bit differences. At least twice that for 8 bit differences.
36	+
37	+	``` pascal
38	+	program JPEGInfo;
39	+
40	+	{$APPTYPE CONSOLE} {$R-}
41	+
42	+	uses strutils,sysutils,classes, md5,magick_wand, ImageMagick,DateUtils;
43	+
44	+	{ Bytes = array[0..0] of Byte;}
45	+	{ LongWords = array[0..255] of Byte;}
46	+
47	+	type
48	+	clr = record
49	+	red:byte;
50	+	green:byte;
51	+	blue:byte;
52	+	opacity:byte;
53	+	end;
54	+	var
55	+	allfile:text;
56	+	infile:text;
57	+	donedir,donefname, inrec:string;
58	+	skipped,deleted,added,modified,count:longword;
59	+	wand: PMagickWand;
60	+	updatefileexists, skipmd5check,notstarted,readnextrec :boolean;
61	+	dctmatrix : array[0..7,0..31] of real;
62	+	status: MagickBooleanType;
63	+	starttime:tdatetime;
64	+	filename:pchar;
65	+	udir,ufname:string;
66	+	Size : Integer;
67	+	timed:longint;
68	+	DirName: string;
69	+
70	+	procedure ThrowWandException(wand: PMagickWand);
71	+	var
72	+	description: PChar;
73	+	severity: ExceptionType;
74	+	begin
75	+	description := MagickGetException(wand, @severity);
76	+	WriteLn(Format('An error occurred. Description: %s', [description]));
77	+	description := MagickRelinquishMemory(description);
78	+	Abort;
79	+	end;
80	+
81	+	function pad(msg:string;len:byte):string;
82	+	var par : byte;
83	+	begin
84	+	result := msg;
85	+	if length(msg) > len then result := copy(msg,1,len);
86	+	while length(result) < len do result := result + ' ';
87	+	for par := 1 to len do if ord(result[par]) < 32 then result[par] := ' ';
88	+	end;
89	+
90	+	{read the next old file record into inrec and fileid into donefname,donedir}
91	+	procedure getnextok;
92	+	var i: integer;
93	+
94	+	procedure skipitem;
95	+	begin
96	+	while (i < length(inrec)) and (inrec[i] <> ' ') do inc(i);
97	+	end;
98	+
99	+	procedure skipspaces;
100	+	begin
101	+	while inrec[i] = ' ' do inc(i)
102	+	end;
103	+	begin
104	+	if eof(infile) then donedir := chr(255)
105	+	else begin
106	+	readln(infile,inrec);
107	+	i := 1;
108	+	skipitem; {md5} skipspaces; {fb54e12c44d2b14d4d3266251109799e}
109	+	skipitem; {md6} skipspaces; {2318170444932902d0541895584df4c2f8f38ad8b1d4dad899870760817422f2}
110	+	skipitem; {phash} skipspaces; {CCFB4EAEF3F30C0D}
111	+	skipitem; {date} skipspaces; {20151001010847}
112	+	skipitem; {size} skipspaces; { 901}
113	+	skipitem; {depth} skipspaces; { 8}
114	+	skipitem; {width} skipspaces; { 32}
115	+	skipitem; {height} skipspaces; { 32}
116	+	skipitem; {name} skipspaces; {01.gif }
117	+	donefname:= copy(inrec,i,999); {F:\CCC\01.gif}
118	+	i := rpos('\',donefname);
119	+	if i = 0 then writeln('error reading line from yalls:',inrec);
120	+	donedir := StringReplace(uppercase(copy(donefname,1,i)),'\',' ',[rfReplaceAll]);;
121	+	donefname := uppercase(copy(donefname,i+1,9999));
122	+	end;
123	+	readnextrec := false;
124	+	end;
125	+	function GetDif(time:longint):string;
126	+	var part: longint; t:longint;
127	+	begin
128	+	result := '';
129	+	t := time;
130	+	if time > 86400 then begin
131	+	part := t div 86400;
132	+	t := t mod 86400;
133	+	result := result + format('%.2d',[part])+' '; end;
134	+	if time > 3600 then begin
135	+	part := t div 3600;
136	+	t := t mod 3600;
137	+	result := result + format('%.2d',[part])+':'; end;
138	+	if time > 60 then begin
139	+	part := t div 60;
140	+	t := t mod 60;
141	+	result := result + format('%.2d',[part])+':'; end;
142	+	if time > 0 then
143	+	result := result + format('%.2d',[t]);
144	+	end;
145	+
146	+	function gethash:string;
147	+	var
148	+	intermediate : array[0..7,0..31] of real;
149	+	i,j,r,c,k:byte;
150	+	img:array[0..31,0..31] of byte;
151	+	dctresult:array[0..7,0..7] of real;
152	+	pimg: Pimage;
153	+	pack:record case integer of
154	+	1:( pkt:PPixelPacket);
155	+	2:( co:^clr);
156	+	end;
157	+	hash,one:qword; mean:DOUBLE;
158	+	begin
159	+	result := 'xxxxxxxxxxxxxxxx'; { in case an error occurs}
160	+	try
161	+	{Resize the image to 32x32. If it is my test case file
162	+	'f:\ccc\01.gif' then it is already at the correct size and I
163	+	found that if imagemagic resizes it the file gets changed ---
164	+	I did not want that, it might be faster to grayscale the image
165	+	first then resize but not if you use MagickQuantizeImage too slow}
166	+	if uppercase(filename) <> 'F:\CCC\01.GIF' then
167	+	MagickResizeImage(wand, 32, 32, LanczosFilter, 1.0);
168	+	{now we convert the image to grayscale - this procedure is similar to
169	+	MagickQuantizeImage(wand, 256, GRAYColorspace, iszero, iszero, iszero);
170	+	but is a LOT faster. BTW I don't know why I had to use variables
171	+	for this function. The conversion numbers were selected to be
172	+	close to what MagickQuantizeImage used. Note I had an issue
173	+	here. Under Q8 IM GetAuthenticPixels it returns an array of
174	+	byte values for colors but the lazarus interface mapped it to word sizes.
175	+	You could get an entire array of values with getauthenticpixels but
176	+	the speed won't be much different so I didn't bother.
177	+	We dump the results into the img array}
178	+	pimg := GetImageFromMagickWand(wand);
179	+	for j := 0 to 31 do begin
180	+	for i := 0 to 31 do begin
181	+	pack.pkt := GetAuthenticPixels(pimg, i, j, 1, 1, nil);
182	+	img[i,j] :=round(pack.co^.red*0.072+
183	+	pack.co^.green*0.715+
184	+	pack.co^.blue*0.213);
185	+	end;
186	+	end;
187	+	{results are in the img array. Now do matrix multiplcation
188	+	DCTMATRIX * img * transposed(DCTMATRIX)
189	+	since we only need the top 8x8 array we only do the
190	+	calculations for that part}
191	+	{intermediate = DCTMATRIX*IMG}
192	+	fillchar(intermediate,sizeof(intermediate),0);
193	+	for r := 0 to 7 do
194	+	for c := 0 to 31 do
195	+	for k := 0 to 31 do
196	+	intermediate[r,c] := intermediate[r,c] +
197	+	dctmatrix[r,k]*img[k,c];
198	+	{dctresult = intermediate*transposed(DCTMATRIX)}
199	+	fillchar(dctresult,sizeof(dctresult),0);
200	+	for r := 0 to 7 do
201	+	for c := 0 to 7 do
202	+	for k := 0 to 31 do
203	+	dctresult[r,c] := dctresult[r,c] +
204	+	intermediate[r,k]*dctmatrix[c,k]{inverted dctmatrix};
205	+	{now we need to calculate the arithmetic mean. Note the first value
206	+	(0,0) is much too large and affects the mean value too much.
207	+	To fix it is set to 0 but perhaps some other solution can be done
208	+	such as do a log on the value}
209	+	mean := 0;
210	+	dctresult[0,0] := 0;
211	+	for r := 0 to 7 do
212	+	for c := 0 to 7 do
213	+	mean := mean + dctresult[r,c];
214	+	mean := mean / 64;
215	+	{the hash is calculated by setting the bit of
216	+	each value greater than the mean}
217	+	one := 1;
218	+	hash := 0;
219	+	for r := 0 to 7 do
220	+	for c := 0 to 7 do begin
221	+	if dctresult[r,c] > mean then inc(hash,one);
222	+	one := one shl 1;
223	+	end;
224	+	one := 1;
225	+	{return the phash}
226	+	result := inttohex(hash,16);
227	+	finally
228	+	end;
229	+	end;
230	+
231	+	procedure processfile(dir:string;fname:string);
232	+	var
233	+	F: file;
234	+	idents:pchar;
235	+	i:integer;
236	+	md5,md6, date, fsize,depth,imageheight,imagewidth,fno:string;
237	+	P: Pointer;
238	+	begin
239	+	inc(count);
240	+	if (count mod 1000) = 1 then begin
241	+	timed := SecondsBetween(Now,starttime);
242	+	writeln('at count ',count-1,' E:'+ GetDif(timed)+' D:',deleted,' A:',added,' M:',modified,' S:',skipped,' ',dir,fname);
243	+	end;
244	+	filename := pchar(dir+fname);
245	+	udir := StringReplace(uppercase(dir),'\',' ',[rfReplaceAll]);;
246	+	ufname := uppercase(fname);
247	+	if updatefileexists then begin
248	+	if readnextrec then getnextok;
249	+	while (udir > donedir) or ((udir = donedir) and (ufname > donefname)) do begin getnextok; inc(deleted); end;
250	+	if skipmd5check then
251	+	if (udir+ufname) = (donedir+donefname) then begin
252	+	inc(skipped);
253	+	writeln(allfile,inrec);
254	+	readnextrec := true;
255	+	exit;
256	+	end;
257	+	end;
258	+	TRY
259	+	Assign(F,FileName);
260	+	FileMode:=fmOpenRead + fmShareDenyNone;
261	+	Reset(F, 1);
262	+	Size := FileSize(F);
263	+	GetMem(P, Size);
264	+	BlockRead(F, P^, Size);
265	+	md5 := md5print(md5buffer(p^,Size));
266	+	finally
267	+	FreeMem(P);
268	+	CloseFile(F);
269	+	end;
270	+	if updatefileexists then begin
271	+	if (udir+ufname) = (donedir+donefname) then begin
272	+	if copy(inrec,1,32) = md5 then begin
273	+	inc(skipped);
274	+	writeln(allfile,inrec);
275	+	readnextrec := true;
276	+	exit;
277	+	end;
278	+	inc(modified);
279	+	end
280	+	else inc(added);
281	+	end;
282	+	if notstarted then begin
283	+	notstarted := false;
284	+	writeln('Started at line:',count,' with file:',filename);
285	+	end;
286	+	try
287	+	wand := NewMagickWand;
288	+	status := MagickReadImage(wand, filename);
289	+	if (status = MagickFalse) then ThrowWandException(wand);
290	+	idents := magickidentifyimage(wand);
291	+	md6 := copy(idents,pos('signature:',idents)+11,64);
292	+	i := pos('date:modify:',idents)+13;
293	+	if i = 0 then i := pos('date:create:',idents)+13;
294	+	date := copy(idents,i ,4)+copy(idents,i+ 5,2)+copy(idents,i+ 8,2)+
295	+	copy(idents,i+11,2)+copy(idents,i+14,2)+copy(idents,i+17,2);
296	+	fsize := copy(idents,pos('Filesize:',idents)+10,64);
297	+	i := POS(Chr(10),FSIZE);
298	+	if i > 0 then fsize := copy(fsize,1,i-1);
299	+	i := pos('Depth:', idents)+7;
300	+	depth := copy(idents, i,1);
301	+	i := pos('Geometry:',idents)+10;
302	+	imageheight := copy(idents, i,99);
303	+	i := pos('x',imageheight);
304	+	imagewidth := copy(imageheight,i+1,99);
305	+	imageheight := copy(imageheight,1,i-1);
306	+	i := pos('+',imagewidth);
307	+	imagewidth := copy(imagewidth,1,i-1);
308	+	fno := fname;
309	+	if length(fno) < 22 then fno := pad(fno,25);
310	+	writeln(allfile,md5,' ',md6,' ',gethash,' ',date:14,' ',
311	+	size:12, depth:2,imageheight:7,imagewidth:7,' ',
312	+	fno,' ',filename);
313	+	finally
314	+	wand := DestroyMagickWand(wand);
315	+	end;
316	+	end;
317	+
318	+	// Recursive procedure to build a list of files
319	+	procedure FindFiles(StartDir: string);
320	+	var
321	+	SR: TSearchRec;
322	+	IsFound: Boolean;
323	+	begin
324	+	if StartDir[length(StartDir)] <> '\' then StartDir := StartDir + '\';
325	+	{ Build a list of the files in directory StartDir
326	+	(but not the directories!) }
327	+	IsFound := FindFirst(StartDir+'.', faAnyFile-faDirectory, SR) = 0;
328	+	while IsFound do begin
329	+	if (uppercase(extractFileExt(sr.name)) = '.JPG')
330	+	or (uppercase(extractFileExt(sr.name)) = '.GIF')
331	+	then processfile(StartDir,SR.Name);
332	+	IsFound := FindNext(SR) = 0;
333	+	end;
334	+	FindClose(SR);
335	+	// Build a list of subdirectories
336	+	IsFound := FindFirst(StartDir+'.', faDirectory, SR) = 0;
337	+	while IsFound do begin
338	+	if ((SR.Attr and faDirectory) <> 0) and (SR.Name[1] <> '.') then
339	+	FindFiles(StartDir + SR.Name);
340	+	IsFound := FindNext(SR) = 0;
341	+	end;
342	+	FindClose(SR);
343	+	end;
344	+
345	+	{$R *.res}
346	+	procedure makedct;
347	+	var c1:real; ii,p,q: byte;
348	+	begin
349	+	c1 := 1/sqrt(32);
350	+	for ii := 0 to 31 do dctmatrix[0,ii] := c1;
351	+	c1 := c1*2;
352	+	for p := 1 to 7 do
353	+	for q := 0 to 31 do
354	+	dctmatrix[p,q] := c1 * cos(pi/64p(2*q+1));
355	+	end;
356	+
357	+	begin
358	+	starttime := now; {keep track of time}
359	+	{have we started processing files as opposed to copy old recs}
360	+	notstarted := true;
361	+	skipmd5check := false; {on copy skip md5 check}
362	+	DirName := ParamStr(1); {did he pass a filename}
363	+	if length(dirname) > 3 then
364	+	if copy(upcase(dirname),1,4) = 'SKIP' then begin
365	+	dirname := paramstr(2);
366	+	skipmd5check := true;
367	+	end;
368	+	count := 0;
369	+	deleted := 0; modified := 0; added := 0; skipped := 0;
370	+	{allow access to file while writing it}
371	+	FileMode:=fmOpenwrite + fmShareDenyNone;
372	+	assign(allfile,'yall');
373	+	rewrite(allfile);
374	+	readnextrec := true;
375	+	updatefileexists := fileexists('yalls');
376	+	if updatefileexists then begin
377	+	{allow access to file while reading it}
378	+	FileMode:=fmOpenwrite + fmShareDenyNone;
379	+	assign(infile,'yalls');
380	+	reset(infile);
381	+	getnextok;
382	+	end;
383	+
384	+	MagickWandGenesis;
385	+	makedct; {make DCT matrix (or at least what is needed)}
386	+	{This is a set of directories that I desired to maintain}
387	+	if (dirname = 'F:\') or (dirname = 'F:\') or
388	+	(length(dirname) < 3) then begin
389	+	FindFiles('F:\CCC\');
390	+	FindFiles('F:\CD\');
391	+	FindFiles('F:\CHECK\');
392	+	FindFiles('F:\QDISK\');
393	+	end
394	+	else FindFiles(DirName); {or use a user supplied directory name}
395	+	writeln('found:',count,' files ',' Deleted:',deleted,' Added:',added,' Changed:',modified,' SKIPPED:',skipped);
396	+	close(allfile);
397	+	if updatefileexists then close(infile);
398	+	MagickWandTerminus;
399	+	{ I get paranoid and like to keep a few versions of the files around}
400	+	if fileexists('yalls.005') then deletefile('yalls.005');
401	+	if fileexists('yalls.004') then renamefile('yalls.004','yalls.005');
402	+	if fileexists('yalls.003') then renamefile('yalls.003','yalls.004');
403	+	if fileexists('yalls.002') then renamefile('yalls.002','yalls.003');
404	+	if fileexists('yalls.001') then renamefile('yalls.001','yalls.002');
405	+	if fileexists('yalls') then renamefile('yalls','yalls.001');
406	+	renamefile('yall','yalls');
407	+	end.
408	+	```
409	+
410	+	---
411	+	From the Free Pascal wiki, page "PHASH", as it stood on 2022-12-11 after a change by Alextp, out of the wiki's dump of 2024-05-26. Copyright stays with the author who wrote it, according to the old wiki's [[Licensing]] page. If you have concerns about this material, [message the admin](/wiki-contact).
412	+	Category: Code
413	+	Category: FPC
414	+	Category: Graphics
415	+	Category: Software security

Back to history