PHASH: what changed

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

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

Back to history