snibgo's ImageMagick pages

Perceptual hash tests

Experiments with IM's perceptual hash, with suggestions for improvement.

lena_rot_-45.pngjpg

This page was developed in response to a forum thread, Get perceptual hash value for image using command line.

References

Phash calculation

Comparing a pair of images by perceptual hash involves two calculations:

  1. Calculate a number of perceptual hash values for each image, individually (based on Hu invariant moments).
  2. Calculate the RMS of the differences between the value sets.

The first calculation is by far the most expensive. The second calculation has variations: a diffferent type of difference, sum instead of mean, not taking the square root, and so on.

In the current implementaton of ImageMagick, seven perceptual hash (PH) values are calculated per channel. There are caclulated in two colorspaces (sRGB and HCL). Each of these colorspaces have three channels, totalling six channels, so there are 42 values. When requested by identify or convert or compare, these 42 values are calculated. When images are compared with -metric phash, an RMS or similar process is applied to the two sets of values.

The process works fairly well. However, for performance, when we are comparing one image to a series of images we don't want to re-calculate the PH values of the same image multiple times, so instead of using -metric phash we should read the PH values directly and then do our own RMS calculation in a script or program.

This raises questions such as:

  1. How many colorspaces are best?
  2. What colorspaces should be used?
  3. Should we use hash values from all the channels?

Objectives

My objectives for perceptual hash comparisons are:

  1. Every pair of images has a score that is zero when the images are identical, with larger values indicating greater perceptual differences.
  2. Images that are essentially the same will have a score that is below a certain threshold.
  3. Images that are not essentially the same will have a score that is above a certain threshold.
  4. Given a new image, we can readily find images that are identical or similar. These will be ordered by smilarity, and the order will be more or less what a human would give.

Variables

The following experiments test the effect of changing:

  1. One colorspace or more.
  2. Choice of colorspace(s).
  3. Use an H?? colorspace but ignore the perceptual hash values from the hue channel.
  4. Use an Y?? colorspace but ignore the perceptual hash values from the Y channel.
  5. Simple difference or proportional diference.

Methodology

We take a small set of standard images, and distort each one in a number of ways, making a dataset 656 images. We use IM to convert each image into 30 different colourspaces and calculate the PH for each. We combine those results to make PH sets for all combinations of two, three or four colorspaces.

Then we calculate RMS differences between PH values for all pairs of images, in each colorspace combination. Comparing error rates tells us the best combinations.

The script setPhEnv.bat sets up the environment:

call %PICTBAT%setPhEnv

Image variations

The standard images are:

dir %DATASET1_DIR% 
 Volume in drive F is New Volume
 Volume Serial Number is F8C5-053E

 Directory of f:\web\im\testImages

09/08/2016  19:57    <DIR>          .
09/08/2016  19:57    <DIR>          ..
09/08/2016  19:55           420,551 airplane.png
09/08/2016  19:56           515,960 barbara.png
09/08/2016  19:56           424,608 boats.png
09/08/2016  19:56           473,054 lena.png
09/08/2016  19:57           621,627 mandril.png
09/08/2016  19:57           501,910 peppers.png
09/08/2016  19:57           437,052 tiffany.png
09/08/2016  19:57           409,896 zelda.png
               8 File(s)      3,804,658 bytes
               2 Dir(s)  54,025,121,792 bytes free

You can download the zip of the images: pht_data1.zip.

For the methods of attack, I have essentially followed Fred's methods. I have cropped each rotated image to the central 50%. This is because photographic images do not usually have triangles of any colour at their edges. I have added vertical and horizontal crops, and hue-shifts. I haven't included the -distort methods. There are 81 attacks of each standard image, making a dataset of (81+1)*8 = 656 images.

The script tweakImages.bat creates the tweaked ("attacked") images. This takes about 3 minutes.

set DATASET1_DIR=\web\im\testImages

call %PICTBAT%tweakImages %DATASET1_DIR% %PH_TWEAKED%

dir %PH_TWEAKED% | tail 
20/08/2016  13:45           111,193 zelda_rot_-90.png
20/08/2016  13:45           158,924 zelda_rot_15.png
20/08/2016  13:45           199,435 zelda_rot_30.png
20/08/2016  13:45           215,533 zelda_rot_45.png
20/08/2016  13:45           123,732 zelda_rot_5.png
20/08/2016  13:45           111,247 zelda_rot_90.png
20/08/2016  13:45           406,560 zelda_trnsp.png
20/08/2016  13:46           406,662 zelda_trnsv.png
             656 File(s)    528,490,419 bytes
               2 Dir(s)  54,025,121,792 bytes free

Colorspaces

We will generate PHs for these colourspaces, which each have three channels:

echo %COLSP3% 
sRGB RGB HCL HSB HSI HSL HWB Lab Luv LCH LCHuv LMS OHTA xyY XYZ YCbCr YCC YDbDr YIQ YUV 

In addition, we use these colorspaces which effectively have only two channels because we ignore channel 0, which is hue or Y:

echo %COLSP_NO_H% 
echo %COLSP_NO_Y% 
CL SB SI SL WB 
CbCr CC DbDr IQ UV 

Those sets are combined, to form this total set of 30 colorspaces:

echo %COLSPACES% 
sRGB RGB HCL HSB HSI HSL HWB Lab Luv LCH LCHuv LMS OHTA xyY XYZ YCbCr YCC YDbDr YIQ YUV CL SB SI SL WB CbCr CC DbDr IQ UV 

Main calculation

The script calcPh.bat calculates the perceptual hashes (PHs) of the images. This uses convert to convert to the colorspace and extract the first column of perceptual hashes, 21 or 14 numbers, writing them to CSV files, one file per colorspace. Processing the 656 tweaked images takes about 2 minutes for each of the 30 colorspaces.

call %PICTBAT%calcPh %PH_TWEAKED%\*

The next processing is more complex. The colourspace combinations for two or more colorspaces are made by simply joining text records from the appropriate single-colorspace CSV files. Then, from the CSV file of PH values for each colorspace combination, the program scorePairs.exe calculates the RMS of differences, comparing every record with every other record. This program has two versions, either with or without proportional difference calculation.

The script ph1colsp.bat writes results for each colorspace and combination of multiple colorspaces to a CSV file, one line per image. It calls the program scorePairsPD.c or scorePairsNPD.c to compare values from all the images pair-wise with each other.

There are 30 single colorspaces, 435 combinations of two colorspaces, 4060 combinations of three colorspaces, and 27405 combinations of four colorspaces, for a total of 31,930 combinations. There is one CSV file per combination, totalling 12 GB of data.

There are 656 variations of images (the originals plus the tweaks), so each colorspace combination requires 656*655/2=214,840 comparisons, for a total of 6,859,841,200 comparisons. This takes a few hours.

(The source data for the comparison, the PH values, are stored as text in CSV files. The scorePairs program re-reads this each time is is run. Storing the data in a binary file would save time.)

call %PICTBAT%ph1colsp %PH_TWEAKED%\*

Results

We can hope that, for each of the 656 images:

  1. The nearest image is in the same group. For example, we hope that the closest image to zelda_rot_45.png is another image in the zelda group. Count the failures: nErrNMLo.
  2. The nearest standard image is the image that was the actual source for this variation. For example, we hope that the closest standard image to zelda_rot_45.png is zelda.png. Count the failures: nErrNrSrc.

The second test is tougher than the first.

These two nearest-neighbour failure rates are shown as the first two columns in the following tables. The most successful colorspace, or combination of colorspaces, is shown first.

Here are the results when only one colorspace is used, for all colorspaces.

copy /y ph1colsp.lis *.csv
cPrefix /iph1colsp.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colsp.csv /t"One colorspace"
call csv2tabNH ph1colsp
One colorspace
nErrNMLo nErrNrSrc highest hiMat loNonMat
4 55 4.90168 3.87167 0.34487 HSB
5 48 3.97985 2.58716 0.309833 YDbDr
5 65 6.60001 3.64642 0.27204 RGB
6 50 3.11678 1.94779 0.297646 OHTA
7 62 3.11844 2.05706 0.294928 YIQ
7 62 5.0583 3.91584 0.383549 HSI
8 55 4.97347 4.49268 0.282611 HCL
8 62 3.18244 2.0729 0.222654 Lab
8 73 3.79216 2.79558 0.241437 LCHuv
8 102 6.00152 3.92534 0.309056 HWB
9 49 3.29632 1.91103 0.31873 Luv
9 66 5.39107 4.83925 0.336949 HSL
11 63 3.16799 1.97398 0.321614 YUV
11 64 3.17103 1.98734 0.306166 YCbCr
12 64 4.31027 2.26602 0.376184 xyY
13 79 3.25221 1.97694 0.239928 YCC
13 100 5.36253 3.40555 0.232492 SI
16 77 3.64582 2.65763 0.226883 LCH
16 82 6.81331 3.56057 0.184021 XYZ
16 82 6.8851 3.61366 0.186627 LMS
17 70 4.99262 3.06772 0.309925 sRGB
17 87 4.92204 3.33453 0.171062 SB
18 62 5.07288 4.34862 0.240966 CL
19 64 4.51441 2.91587 0.218918 DbDr
23 108 5.59984 4.87462 0.133584 SL
24 77 3.26735 2.2019 0.187485 IQ
28 137 7.28995 4.79995 0.217739 WB
32 60 2.79793 2.05945 0.176273 UV
33 67 2.73093 2.11339 0.14736 CC
33 67 2.8137 2.07207 0.146332 CbCr

The table also shows the highest (worst) score achieved by one image against another. The best score is always 0.0, when two images are identical.

hiMat is the highest (worst) score between two images that should match, eg two versions of Zelda. loNonMat is the lowest score between two images that shouldn't match, eg a version of Zelda and a version of Barbara. Ideally, hiMat would be less than loNonMat. Then we could pick any number between these two as a threshold for deciding whether images were different. Sadly this isn't true, we have a negative margin, so any number we pick as a threshold will find either false positives (incorrectly declaring two images are the same) or false negatives (incorrectly declaring two images are different) or both.

The attacks lead to a wide variation in versions of an image, and hiMat gives the worst score between extremes.

Are two colorspaces better than one? Of the 435 results, we show the best 15.

cHead /iph1colsp2.lis /oph1colsp2.csv /h15
cPrefix /iph1colsp2.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colsp2.csv /t"Two colorspaces"
call csv2tabNH ph1colsp2
Two colorspaces
nErrNMLo nErrNrSrc highest hiMat loNonMat
0 20 3.66768 2.7596 0.608886 HSB xyY
0 23 3.77649 2.22186 0.484251 xyY YDbDr
0 25 3.53391 2.73898 0.453201 HSB YIQ
0 25 3.5666 2.74955 0.513635 HSB YUV
0 26 3.89057 3.48117 0.425032 HCL IQ
0 27 3.76144 2.3041 0.47228 DbDr xyY
0 27 5.11442 2.85463 0.51817 HSB RGB
0 28 3.51981 2.74309 0.54539 HSB Lab
0 28 3.56357 2.74999 0.513971 HSB YCbCr
0 28 3.56499 2.74625 0.555173 HSB Luv
0 28 3.59668 3.17792 0.39523 HCL YIQ
0 28 3.77271 2.75139 0.514552 HSB YDbDr
0 28 3.99317 3.1957 0.52362 HCL xyY
0 28 4.16309 2.73912 0.526136 HSB sRGB
0 29 3.94899 3.43942 0.535023 HSL xyY

IM currently uses colorspaces HCL and sRGB. The score for this pair of colorspaces is:

nErrNMLo nErrNrSrc highest hiMat loNonMat
3 32 4.10509 3.17804 0.296583 HCL sRGB

This pair of colorspaces performs better than either alone, but worse than many other pairs.

Are three colorspaces better than two? Of the 4060 results, we show the best 15.

cHead /iph1colsp3.lis /oph1colsp3.csv /h15
cPrefix /iph1colsp3.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colsp3.csv /t"Three colorspaces"
call csv2tabNH ph1colsp3
Three colorspaces
nErrNMLo nErrNrSrc highest hiMat loNonMat
0 12 3.41657 2.39093 0.588125 HSB IQ xyY
0 13 3.35277 2.39991 0.592992 CC HSB xyY
0 13 3.35608 2.40039 0.615952 CbCr HSB xyY
0 13 3.35914 2.40001 0.621343 HSB UV xyY
0 13 3.45529 2.29031 0.59955 HSB xyY YDbDr
0 13 4.53723 2.57957 0.638959 HWB xyY YDbDr
0 14 3.14664 2.3037 0.55124 HSB YDbDr YIQ
0 14 3.36009 2.4016 0.61108 DbDr HSB xyY
0 14 3.40652 2.41775 0.611731 HSI IQ xyY
0 14 3.48114 2.25425 0.578785 HSB xyY YIQ
0 14 3.48226 2.26014 0.627531 HSB Luv xyY
0 14 3.48838 2.27954 0.600141 HSI xyY YIQ
0 15 3.45182 2.25434 0.602629 HSB OHTA xyY
0 15 3.45293 2.26272 0.583164 HSB xyY YCC
0 15 3.46084 2.25757 0.579672 HSB Lab xyY

Are four colorspaces better than three? Of the 27,405 results, we show the best 15.

cHead /iph1colsp4.lis /oph1colsp4.csv /h15
cPrefix /iph1colsp4.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colsp4.csv /t"Four colorspaces"
call csv2tabNH ph1colsp4
Four colorspaces
nErrNMLo nErrNrSrc highest hiMat loNonMat
0 9 4.19069 2.29402 0.605738 HWB xyY YDbDr YIQ
0 10 3.28461 2.26171 0.582039 DbDr HSB xyY YIQ
0 10 3.28461 2.26171 0.582039 HSB IQ xyY YDbDr
0 10 3.32219 2.14946 0.600949 HSB xyY YDbDr YUV
0 10 3.33027 2.17089 0.576159 HSB Lab xyY YDbDr
0 10 3.34611 2.2009 0.57549 HSB xyY YDbDr YIQ
0 10 3.34698 2.15893 0.638722 HSB Luv xyY YDbDr
0 10 4.07347 2.4541 0.628747 DbDr HWB IQ xyY
0 10 4.45469 2.41191 0.647211 HSB RGB xyY YDbDr
0 10 4.65953 2.58858 0.558847 DbDr HSB LMS sRGB
0 10 4.67569 2.57206 0.542437 HSB LMS sRGB YDbDr
0 11 3.19735 2.1718 0.612624 CbCr HSI IQ xyY
0 11 3.2562 2.07431 0.605115 CbCr HSB OHTA xyY
0 11 3.27702 2.24538 0.599474 DbDr HSI xyY YIQ
0 11 3.27702 2.24538 0.599474 HSI IQ xyY YDbDr

In the lists above, what colorspaces occur most frequently?

call %PICTBAT%phUnion pht_popcolsp.csv
call csv2tabNH pht_popcolsp
Colorspace Count
xyY 33
HSB 32
YDbDr 15
YIQ 10
IQ 8
DbDr 7
HSI 6
sRGB 4
Luv 4
CbCr 4
HCL 4
HWB 4
Lab 4
OHTA 3
YUV 3
LMS 3
RGB 3
YCbCr 2
HSL 2
YCC 2
UV 2
CC 2
XYZ 1
WB 1
SL 1
SI 1
SB 1
LCHuv 1
LCH 1
CL 1

Experiments with different data suggest that this list is more stable than the order of colorspaces in the individual tests.

The set of xyY, HSB, YDbDr and YIQ seem a suitable set.

Frankly, I'm surprised that HSB performs so well. In other tests where a greater proportion of attacks have hue shifts, it performs less well.

Proportional difference

The usual RMS of differences uses simple differences. That is, it assumes that the difference between 10.0 and 10.5 has the same significance as the difference between 1.0 and 1.5. Perhaps this is wrong, and it is as significant as the difference between 1.00 and 1.05.

Dividing the difference by the sum (or average) of the values is one way to compensate. However, values can be negative, so the sum can be close to or exactly zero, over-exagerating the difference. I set the near-zero bound of the sum at 0.25. This value was found empirically to minimise errors.

This is implemented in scorePairs.inc, in the section #if PROP_DIFF==1 ... #endif.

To show the improvement from proportional differences, we run the test comparing 656 images with each other on all colorspaces, first using the proportional difference, then the simple difference.

call %PICTBAT%ph1colspPD %PH_TWEAKED%\*
if ERRORLEVEL 1 goto :error

call %PICTBAT%ph1colspNPD %PH_TWEAKED%\*
cSort /iph1colspPD.lis /oph1colspPD.csv /k5
cPrefix /iph1colspPD.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colspPD.csv /t"Proportional difference"
call csv2tabNH ph1colspPD

cSort /iph1colspNPD.lis /oph1colspNPD.csv /k5
cPrefix /iph1colspNPD.csv /tnErrNMLo,nErrNrSrc,highest,hiMat,loNonMat
cPrefix /iph1colspNPD.csv /t"Not PD"
call csv2tabNH ph1colspNPD

These results are sorted by the colorspace, for direct comparison between the methods.

Proportional difference
nErrNMLo nErrNrSrc highest hiMat loNonMat
31 65 0.258922 0.21723 0.0170672 CC
8 76 4.83968 4.55564 0.033006 CL
11 41 0.319059 0.228104 0.0252825 CbCr
3 38 1.62197 1.62197 0.0407106 DbDr
5 42 3.97132 3.73921 0.0706143 HCL
3 45 1.27138 1.27138 0.0532312 HSB
3 54 4.49538 4.15535 0.0640453 HSI
5 54 1.43523 1.43523 0.0605576 HSL
4 87 12.3469 11.5198 0.0851665 HWB
13 45 0.306826 0.221677 0.0235304 IQ
4 61 0.661678 0.462428 0.0367203 LCH
6 60 0.736318 0.502837 0.0418211 LCHuv
17 70 2.4591 2.4591 0.0445579 LMS
9 67 0.513136 0.358215 0.0276584 Lab
2 48 0.514812 0.358575 0.0360022 Luv
0 50 0.571311 0.393063 0.0347379 OHTA
9 61 8.1569 8.1569 0.0667914 RGB
16 108 1.47922 1.47922 0.0134013 SB
23 121 5.48576 5.06932 0.0341505 SI
16 103 1.68414 1.68414 0.0242487 SL
13 41 0.306025 0.221876 0.0253182 UV
25 121 15.1085 14.1085 0.0501193 WB
18 78 2.80699 2.80699 0.0435658 XYZ
11 99 0.555052 0.378293 0.0242488 YCC
7 60 0.550152 0.375209 0.0276769 YCbCr
0 33 1.32462 1.32462 0.0580899 YDbDr
3 65 0.551073 0.373374 0.031116 YIQ
7 65 0.550068 0.374502 0.0298989 YUV
0 61 1.0526 0.761962 0.0562753 sRGB
2 54 0.982209 0.685397 0.0541901 xyY
Not PD
nErrNMLo nErrNrSrc highest hiMat loNonMat
33 67 2.73093 2.11339 0.14736 CC
18 62 5.07288 4.34862 0.240966 CL
33 67 2.8137 2.07207 0.146332 CbCr
19 64 4.51441 2.91587 0.218918 DbDr
8 55 4.97347 4.49268 0.282611 HCL
4 55 4.90168 3.87167 0.34487 HSB
7 62 5.0583 3.91584 0.383549 HSI
9 66 5.39107 4.83925 0.336949 HSL
8 102 6.00152 3.92534 0.309056 HWB
24 77 3.26735 2.2019 0.187485 IQ
16 77 3.64582 2.65763 0.226883 LCH
8 73 3.79216 2.79558 0.241437 LCHuv
16 82 6.8851 3.61366 0.186627 LMS
8 62 3.18244 2.0729 0.222654 Lab
9 49 3.29632 1.91103 0.31873 Luv
6 50 3.11678 1.94779 0.297646 OHTA
5 65 6.60001 3.64642 0.27204 RGB
17 87 4.92204 3.33453 0.171062 SB
13 100 5.36253 3.40555 0.232492 SI
23 108 5.59984 4.87462 0.133584 SL
32 60 2.79793 2.05945 0.176273 UV
28 137 7.28995 4.79995 0.217739 WB
16 82 6.81331 3.56057 0.184021 XYZ
13 79 3.25221 1.97694 0.239928 YCC
11 64 3.17103 1.98734 0.306166 YCbCr
5 48 3.97985 2.58716 0.309833 YDbDr
7 62 3.11844 2.05706 0.294928 YIQ
11 63 3.16799 1.97398 0.321614 YUV
17 70 4.99262 3.06772 0.309925 sRGB
12 64 4.31027 2.26602 0.376184 xyY

For almost all colorspaces, the results using the proportional calculation are better.

In addition, the proportional calculation has three colorspaces with a perfect nErrNMLo=0. For the non-proportional, the best score is nErrNMLo=4.

Conclusion: when only one colorspace is used, the proportional difference is better.

However, using the proportional calculation for tests of two, three and four colorspaces has a different result: we again get perfect zeros for nErrNMLo, but the scores for nErrNrSrc are higher (worse):

For two colorspaces, non-proportional has the best 15 nErrNrSrc ranging from 20 to 29, where proportional ranges from 27 to 31.

For three colorspaces, non-proportional has the best 15 nErrNrSrc ranging from 12 to 15, where proportional ranges from 21 to 25.

For four colorspaces, non-proportional has the best 15 nErrNrSrc ranging from 9 to 11, where proportional ranges from 17 to 19.

Conclusion: for phash using multiple colorspaces, the proportional difference is worse.

Ranges of values

IM caps PH values at a maximum of +11.0. In this dataset, about 1% of the values are negative, with no values less than -4.5. There appears to be no clamping of negative values, and values of -100 etc are theoretically possible.

If we take the likely range as -5 to +11, the ordinary difference has a possible absolute range of 0.0 to +16.0.

My proportional difference has a maximum where the inputs are (-5,+5). For those inputs, the difference is (5+5)/0.25 = 40. However, most values will be in the range 0.0 to 1.0.

As I then take the RMS of differences, the possble and likely ranges of the RMS are the same.

Conclusions

No combination of colorspaces will be ideal for every application. There is a clear time versus quality tradeoff: more colorspaces give more accurate results but take longer to calculate.

Two colorspaces are very much better than one. Three or four colorspaces give lesser improvements.

When only one colorspace is used, the choice of colorspace has a large impact on the correctness of the results. (And the proportional difference calculation is better.)

When four colorspaces are used, the choice of colorspace has smaller impact on the correctness of the results, so the colorspace combination is less critical.

The big improvement of two colorspaces instead of one suggests that we should always use at least two.

When two colorspaces are used, they should be xyY and HSB.

When three colorspaces are used, they should be xyY, HSB and IQ. However xyY, HSB and YIQ are nearly as good. (The difference is probably not significant.)

When four colorspaces are used, they should be xyY, HWB, YIQ and YDbDr. However xyY, HSB, YIQ and YDbDr are nearly as good. (Again, the difference is probably not significant.)

I think it likely that future work will refine perceptual hash techniques, and ImageMagick should be engineered to enable this.

Recommendations for ImageMagick internals

For the calculation of perceptual hash values of images, I recommend that:

  1. IM's default behaviour is to calculate PHs for two colorspaces, as at present, but that the default colorspaces should be xyY and HSB.
  2. IM's calculation for differences stays as it is, a simple subtraction.
  3. A facility is provided to use different colorspaces instead of xyY and HSB.
  4. A facility is provided to use any number of colorspaces.
  5. A facility is provided to omit any channels from the colorspaces.

For the calculation of -metric phash comparisons between images, I recommend that:

  1. IM squares the differences and sums them, as at present.
  2. IM then divides by the appropriate number of channels.
  3. IM then takes the square root. Thus we have the conventional RMS calculation.

(For comparative purposes, taking the square root is not important. But taking the correct mean is important.)

IM V7 has (or will have) the capability to add channels to an image. I suppose a command can (or will be able to) add channels for HSB, xyY and so on. So an obvious technique is for IM to calculate perceptual hashes for only the channels in the image, possibly reduced by a -channel option. The application can add as many channels as are required. This would provide all the facilities recommended above.

If that solution isn't feasible, an alternative is to tell IM what colorspaces to use, perhaps with a list, and another option for channel selection:

-define phash:colorspaces=xyY,HSB,YDbDr,YIQ
-define phash:omit_channels=3,6

Instead of selecting channels, we might give abbreviated names in the colorspace list:

-define phash:colorspaces=xyY,SB,DbDr,YIQ

However, this would create ambiguity, eg the U and V channels of Luv and YUV are different.

Future

As a check, we should include tests of two colorspaces that are the same, eg sRGB, and sRGB. We expect the results to be the same as just the one colorspace.

For a given colorspace combination, what are the error rates at various thresholds of RMS? This is addressed in Perceptual hash thresholds.

Experiment with http://www.cis.temple.edu/~latecki/TestData/mpeg7shapeB.tar.gz "Shape data for the MPEG-7 core experiment CE-Shape-1".

Experiment with a dataset of 20,000 images, each up to about 35 M-pixels.

Scripts

For convenience, .bat scripts are also available in a single zip file. See Zipped BAT files.

setPhEnv.bat

rem Set environment for perceptual hash tests.

rem ---------------------------------------
rem Directories. Names must NOT contain underscore.

set DATASET1_DIR=\web\im\testImages

set PH_TWEAKED=phtweaked
if not exist %PH_TWEAKED% md %PH_TWEAKED%

set PH_CSVDIR=phcsv
if not exist %PH_CSVDIR% md %PH_CSVDIR%


rem ---------------------------------------
rem Colorspaces

rem Note:
rem   CL and CLp would give the same result.
rem   YCbCr and YPbPr would give the same result.
rem   RGB and scRGB would give the same result.

set COLSP3=sRGB RGB HCL HSB HSI HSL HWB Lab Luv LCH LCHuv LMS OHTA xyY XYZ YCbCr YCC YDbDr YIQ YUV

rem Next omits channel 0, which is "H".
rem Note:
rem   SB and SV would give the same result.
rem   CL and CLp would give the same result.

set COLSP_NO_H=CL SB SI SL WB

set COLSP_NO_Y=CbCr CC DbDr IQ UV

set COLSPACES=%COLSP3% %COLSP_NO_H% %COLSP_NO_Y%

rem Use all colorspaces:
set COLSP_FOR3=%COLSPACES%

tweakImages.bat

rem Given %1 is a directory containing just test images,
rem %2 is a destination directory,
rem creates tweaked versions of images in %2.

set OUTDIR=%2

set EXT=png

for %%F in (%1\*) do (
  echo %%F

  set BASE=%%~nF

  copy /y %%F %OUTDIR%\%%~nxF

  for %%V in (-20 -15 -10 -5 5 10 15 20) do %IM%convert ^
    %%F ^
    +depth ^
    ^( +clone ^
       -brightness-contrast %%V,0 ^
       +write %OUTDIR%\!BASE!_bri_%%V.%EXT% ^
       +delete ^
    ^) ^
    -brightness-contrast 0,%%V ^
    %OUTDIR%\!BASE!_con_%%V.%EXT%

  for %%V in (0.75 0.8 0.85 0.9 0.95 1.05 1.1 1.15 1.2 1.25) do %IM%convert ^
    %%F ^
    -gamma %%V ^
    -set gamma 0.454545 ^
    +depth ^
    %OUTDIR%\!BASE!_gam_%%V.%EXT%

  for %%V in (0.5 1 2 3) do %IM%convert ^
    %%F ^
    -blur 0x%%V ^
    +depth ^
    %OUTDIR%\!BASE!_blr_%%V.%EXT%

  for %%V in (5 10 25 50 75 100) do %IM%convert ^
    %%F ^
    -quality %%V ^
    %OUTDIR%\!BASE!_jpg_%%V.jpg

  rem Colour gaussian noise.
  for %%V in (0.5 1 1.5 2 2.5 3) do %IM%convert ^
    -seed 1234 ^
    %%F ^
    -attenuate %%V +noise gaussian ^
    +depth ^
    %OUTDIR%\!BASE!_nsc_%%V.%EXT%

  rem Gray gaussian noise.
  for %%V in (0.5 1 1.5 2 2.5 3) do %IM%convert ^
    -seed 1234 ^
    %%F ^
    ^( +clone ^
       -fill gray^(50%%^) -colorize 100 ^
       -attenuate %%V +noise gaussian ^
       -modulate 100,0,100 ^
    ^) ^
    -compose hardlight -composite ^
    +depth ^
    %OUTDIR%\!BASE!_nsg_%%V.%EXT%

  for %%V in (25 50 75 150 200 400) do %IM%convert ^
    %%F ^
    -resize %%V%% ^
    %OUTDIR%\!BASE!_res_%%V.%EXT%

  for %%V in (5 15 30 45 90) do %IM%convert ^
    %%F ^
    -background black ^
    ^( +clone ^
       -rotate %%V +repage ^
       -gravity Center ^
       -crop 50%%x50%%+0+0 +repage ^
       +write %OUTDIR%\!BASE!_rot_%%V.%EXT% ^
       +delete ^
    ^) ^
    -rotate -%%V +repage ^
    -gravity Center ^
    -crop 50%%x50%%+0+0 +repage ^
    %OUTDIR%\!BASE!_rot_-%%V.%EXT%

  for %%V in (90 95 97.5 102.5 105 110) do %IM%convert ^
    %%F ^
    -modulate 100,100,%%V ^
    %OUTDIR%\!BASE!_hue_%%V.%EXT%

  for %%V in (70 80 90) do %IM%convert ^
    %%F ^
    -gravity Center ^
    ^( +clone ^
       -crop %%V%%x100%%+0+0 +repage ^
       +write %OUTDIR%\!BASE!_crpv_%%V.%EXT% ^
       +delete ^
    ^) ^
    -crop 100%%x%%V%%+0+0 +repage ^
    %OUTDIR%\!BASE!_crph_%%V.%EXT%

  %IM%convert ^
    %%F ^
    ^( +clone ^
       -flip ^
       +write %OUTDIR%\!BASE!_flip.%EXT% ^
       +delete ^
    ^) ^
    ^( +clone ^
       -flop ^
       +write %OUTDIR%\!BASE!_flop.%EXT% ^
       +delete ^
    ^) ^
    ^( +clone ^
       -transpose ^
       +write %OUTDIR%\!BASE!_trnsp.%EXT% ^
       +delete ^
    ^) ^
    ^( +clone ^
       -transverse ^
       +write %OUTDIR%\!BASE!_trnsv.%EXT% ^
       +delete ^
    ^) ^
    ^( +clone ^
       -gravity Center -crop 50%%x50%%+0+0 +repage ^
       +write %OUTDIR%\!BASE!-half.%EXT% ^
       +delete ^
    ^) ^
    NULL:

)

rem Note: "half" has hyphen, not underscore.

calcPh.bat

rem rem Given %1 is wildcard file (eg mydir\*.png),
rem writes filename and 21 (or 14) perceptual hash values
rem (first column only)
rem to CSV files for various colorspaces.

set IMAGES=%1

if "%IMAGES%"=="" (
  echo %0: No images.
  exit /B 1
)

call %PICTBAT%setPhEnv

if "%PH_CSVDIR%"=="" exit /B 1

echo %0: %IMAGES%

goto skip

rem Calculate PHs for ordinary colorspaces.
rem

for %%C in (%COLSP3%) do (
  call %PICTBAT%idPh21imgs %IMAGES% %PH_CSVDIR%\idPh21_%%C.csv %%C

  if ERRORLEVEL 1 exit /B 1
)

rem Calculate PHs for H** colorspaces, without the H channel.
rem

for %%C in (%COLSP_NO_H%) do (
  call %PICTBAT%idPh21imgs %IMAGES% %PH_CSVDIR%\idPh21_%%C.csv H%%C 12

  if ERRORLEVEL 1 exit /B 1
)

:skip

rem Calculate PHs for Y** colorspaces, without the Y channel.
rem

for %%C in (%COLSP_NO_Y%) do (
  call %PICTBAT%idPh21imgs %IMAGES% %PH_CSVDIR%\idPh21_%%C.csv Y%%C 12

  if ERRORLEVEL 1 exit /B 1
)

exit /B 0

ph1colsp.bat

rem rem Given %1 is wildcard file (eg mydir\*.png),
rem writes filename and 21 (or 14) perceptual hash values
rem (first column only)
rem to CSV files for various colorspaces.

set IMAGES=%1

if "%IMAGES%"=="" (
  echo %0: No images.
  exit /B 1
)


call %PICTBAT%setPhEnv

if "%PH_CSVDIR%"=="" exit /B 1

set SCORE_PAIRS=scorePairsNPD.exe


rem --------------------------------------------------------------
rem
rem Build the CSV files (idPh21_*.csv) for single colorspaces.

rem for %%C in (%COLSPACES%) do (
rem   call %PICTBAT%idPh21imgs %IMAGES% %PH_CSVDIR%\idPh21_%%C.csv %%C
rem 
rem   if ERRORLEVEL 1 exit /B 1
rem )

set LIST1=ph1colsp.lis

del %LIST1%

for %%C in (%COLSPACES%) do (
  %SCORE_PAIRS% %PH_CSVDIR%\idPh21_%%C.csv | cgrep /p0 /i- /snErr /f"\r,%%C\n" /O%LIST1%

  rem %SCORE_PAIRS% %PH_CSVDIR%\idPh21_%%C.csv | cgrep /p0 /i- /n[ /s"[NonLo[**[" /f"%%C \s\r\n" /O%LIST2%
)

chStrs /i%LIST1% /s%PICTBAT%\ph1colsp.dat
chStrs /i%LIST1% /f\q
csort  /i%LIST1% /o%LIST1% /k0-2

type %LIST1%

rem --------------------------------------------------------------
rem
rem Build the CSV files (idPh42_*_*.csv) for two colorspaces.

set LIST2=ph1colsp2.lis
del %LIST2%

rem Note: cJoin here makes an output from two inputs, appending the numeric columns.
rem (It also sorts the inputs.)

for %%C in (%COLSPACES%) do (
  set INF0=%PH_CSVDIR%\idPh21_%%C.csv
  for %%D in (%COLSPACES%) do (
    if /I %%C LSS %%D (
      set INF1=%PH_CSVDIR%\idPh21_%%D.csv
      set OUTCSV=%PH_CSVDIR%\idPh42_%%C_%%D.csv
      cJoin /p0 /i!INF0!,!INF1! /k0 /K0 /n /o!OUTCSV!

      %SCORE_PAIRS% !OUTCSV! | cgrep /p0 /i- /snErr /f"\r,%%C,%%D\n" /O%LIST2%

    )
  )
)

chStrs /i%LIST2% /s%PICTBAT%\ph1colsp.dat
chStrs /i%LIST2% /f\q
cSort  /i%LIST2% /o%LIST2% /k0-3



rem --------------------------------------------------------------
rem
rem Build the CSV files (idPh63_*_*_*.csv) for three colorspaces.

set LIST3=ph1colsp3.lis
del %LIST3%

set TMPCSV=ph1_tmp.csv

echo off
for %%C in (%COLSP_FOR3%) do (
  set INF0=%PH_CSVDIR%\idPh21_%%C.csv
  for %%D in (%COLSP_FOR3%) do (
    if /I %%C LSS %%D (
      set INF1=%PH_CSVDIR%\idPh21_%%D.csv
      for %%E in (%COLSP_FOR3%) do (
        if /I %%D LSS %%E (
          set INF2=%PH_CSVDIR%\idPh21_%%E.csv
          echo %%C %%D %%E

          set OUTCSV=%PH_CSVDIR%\idPh63_%%C_%%D_%%E.csv
          cJoin /p0 /i!INF0!,!INF1! /k0 /K0 /n /o%TMPCSV%
          cJoin /p0 /i%TMPCSV%,!INF2! /k0 /K0 /n /o!OUTCSV!

          %SCORE_PAIRS% !OUTCSV! | cgrep /p0 /i- /snErr /f"\r,%%C,%%D,%%E\n" /O%LIST3%
        )
      )
    )
  )
)
echo on

del %TMPCSV%

chStrs /i%LIST3% /s%PICTBAT%\ph1colsp.dat
chStrs /i%LIST3% /f\q
cSort  /i%LIST3% /o%LIST3% /k0-4

rem --------------------------------------------------------------
rem
rem Build the CSV files (idPh84_*_*_*_*.csv) for four colorspaces.

set LIST4=ph1colsp4.lis
del %LIST4%

set TMPCSV=ph1_tmp.csv

echo off
for %%C in (%COLSP_FOR3%) do (
  set INF0=%PH_CSVDIR%\idPh21_%%C.csv
  for %%D in (%COLSP_FOR3%) do (
    if /I %%C LSS %%D (
      set INF1=%PH_CSVDIR%\idPh21_%%D.csv
      for %%E in (%COLSP_FOR3%) do (
        if /I %%D LSS %%E (
          set INF2=%PH_CSVDIR%\idPh21_%%E.csv

          for %%F in (%COLSP_FOR3%) do (
            if /I %%E LSS %%F (
              set INF3=%PH_CSVDIR%\idPh21_%%F.csv
              echo %%C %%D %%E %%F

              set OUTCSV=%PH_CSVDIR%\idPh84_%%C_%%D_%%E_%%F.csv
              cJoin /p0 /i!INF0!,!INF1! /k0 /K0 /n /o%OUTCSV%
              cJoin /p0 /i%OUTCSV%,!INF2! /k0 /K0 /n /o!TMPCSV!
              cJoin /p0 /i%TMPCSV%,!INF3! /k0 /K0 /n /o!OUTCSV!

              %SCORE_PAIRS% !OUTCSV! | cgrep /p0 /i- /snErr /f"\r,%%C,%%D,%%E,%%F\n" /O%LIST4%
            )
          )
        )
      )
    )
  )
)
echo on

del %TMPCSV%

chStrs /i%LIST4% /s%PICTBAT%\ph1colsp.dat
chStrs /i%LIST4% /f\q
cSort  /i%LIST4% /o%LIST4% /k0-5



rem --------------------------------------------------------------
rem
rem Build the CSV files (idPh_multi.csv) for many colorspaces.

set LISTMULTI=ph1colspm.lis
del %LISTMULTI%

set TMPCSV=ph1_tmp.csv

set OUTCSV=%PH_CSVDIR%\idPh_multi.csv

copy /Y %PH_CSVDIR%\idPh21_OHTA.csv %OUTCSV%

for %%C in (%COLSP_FOR3%) do (
  if not %%C==OHTA (
    set INF0=%PH_CSVDIR%\idPh21_%%C.csv
    cJoin /p0 /i%OUTCSV%,!INF0! /k0 /K0 /n /o%TMPCSV%
    copy /Y %TMPCSV% %OUTCSV%
  )
)

%SCORE_PAIRS% !OUTCSV! | cgrep /p0 /i- /snErr /f"\r,%COLSP_FOR3%\n" /O%LISTMULTI%


del %TMPCSV%

chStrs /i%LISTMULTI% /s%PICTBAT%\ph1colsp.dat
chStrs /i%LISTMULTI% /f\q
cSort  /i%LISTMULTI% /o%LISTMULTI% /k0-2

exit /B 0

idPh21imgs.bat

rem Given %1 is wildcard file (eg mydir\*.png),
rem writes filename and 21 or 14 or 7 perceptual hash values
rem (first column only)
rem to %2, output CSV file.
rem %3 is colorspace for hashing. (Should have 3 channels.)
rem %4 is which channel to output.
rem   Any of 0, 1, 2. Eg 012 or 12 or 01. Default 012.

setlocal enabledelayedexpansion

set IMAGES=%1

if "%IMAGES%"=="" (
  echo %0: No images.
  exit /B 1
)

set OUTFILE=%2
if "%OUTFILE%"=="." set OUTFILE=
if "%OUTFILE%"=="" set OUTFILE=ip42i.csv

set COLSP=%3
if "%COLSP%"=="." set COLSP=
if "%COLSP%"=="" set COLSP=sRGB

set OUTCH=%4
if "%OUTCH%"=="." set OUTCH=
if "%OUTCH%"=="" set OUTCH=012

set OUT_0=0
set OUT_1=0
set OUT_2=0

set TMP=%OUTCH:0=X%
if not %TMP%==%OUTCH% set OUT_0=1
set TMP=%OUTCH:1=X%
if not %TMP%==%OUTCH% set OUT_1=1
set TMP=%OUTCH:2=X%
if not %TMP%==%OUTCH% set OUT_2=1

if %OUT_0%==0 if %OUT_1%==0 if %OUT_2%==0 (
  echo %0: No channels to output.
  exit /B 1
)

echo %0: %IMAGES% %OUTFILE% %COLSP% %OUT_0% %OUT_1% %OUT_2%

if /I %COLSP%==sRGB (
  set sCOLSP=
) else (
  set sCOLSP=-colorspace %COLSP% -set colorspace sRGB
)


del %OUTFILE% 2>nul

echo off

echo %0: Identify perceptual hashes of files %1 colorspace %COLSP% ...

set nFILES=0
(
  for %%F in (%1) do (
    echo %%F
    set nC=0
    set PH[0][1]B=
    for /F "usebackq tokens=1-3 delims=, " %%A in (`%IM%convert ^
      %%F ^
      -resize "600x600>" ^
      %sCOLSP% ^
      -verbose -define identify:moments ^
      info:`) do (
      rem echo %%A %%B %%C
      if "%%A"=="PH1:" (
        set PH[!nC!][1]B=%%B
        set PH[!nC!][1]C=%%C
      ) else if "%%A"=="PH2:" (
        set PH[!nC!][2]B=%%B
        set PH[!nC!][2]C=%%C
      ) else if "%%A"=="PH3:" (
        set PH[!nC!][3]B=%%B
        set PH[!nC!][3]C=%%C
      ) else if "%%A"=="PH4:" (
        set PH[!nC!][4]B=%%B
        set PH[!nC!][4]C=%%C
      ) else if "%%A"=="PH5:" (
        set PH[!nC!][5]B=%%B
        set PH[!nC!][5]C=%%C
      ) else if "%%A"=="PH6:" (
        set PH[!nC!][6]B=%%B
        set PH[!nC!][6]C=%%C
      ) else if "%%A"=="PH7:" (
        set PH[!nC!][7]B=%%B
        set PH[!nC!][7]C=%%C
        set /A nC+=1
      )
    )

    if "!PH[0][1]B!"=="" exit /B 1

    set sVALS=

    if %OUT_0%==1 set sVALS=!sVALS!^
!PH[0][1]B!,!PH[0][2]B!,!PH[0][3]B!,!PH[0][4]B!,!PH[0][5]B!,!PH[0][6]B!,!PH[0][7]B!,

    if %OUT_1%==1 set sVALS=!sVALS!^
!PH[1][1]B!,!PH[1][2]B!,!PH[1][3]B!,!PH[1][4]B!,!PH[1][5]B!,!PH[1][6]B!,!PH[1][7]B!,

    if %OUT_2%==1 set sVALS=!sVALS!^
!PH[2][1]B!,!PH[2][2]B!,!PH[2][3]B!,!PH[2][4]B!,!PH[2][5]B!,!PH[2][6]B!,!PH[2][7]B!,

    set sVALS=!sVALS:~0,-1!

    echo %%F,!sVALS! >>%OUTFILE%

    set /A nFILES+=1
  )
)

echo nFILES=%nFILES%

endlocal

csv2tabNH.bat

rem From %1.csv, creates partial %1.htm with a table that has no th tags.
rem %2 is parameters to table, eg "border".

set c2tTAB_PARAMS=%2
if "%c2tTAB_PARAMS%"=="." set c2tTAB_PARAMS=

chSep   /p0 /i%1.csv /o%1.htm /f"," /t\(/td\)\(td\) /n
cPrefix /p0 /i%1.htm /l\(tr\)\(td\) /r\(/td\)\(/tr\) /t"\(table %c2tTAB_PARAMS%\)" /b\(/table\) /X

Programs

The scripts also use general-purpose programs such as cPrefix.exe, and I do not supply the source or binaries of these.

scorePairsPD.c

// Whether to use proportional difference.
#define PROP_DIFF 1

#include "scorePairs.inc"

int main (int argc, char *argv [])
{
  BOOL okay = ScorePairs (argv[1], FALSE);

  return (okay ? 0 : 1) ;
}

scorePairsNPD.c

// Whether to use proportional difference.
#define PROP_DIFF 0

#include "scorePairs.inc"

int main (int argc, char *argv [])
{
  BOOL okay = ScorePairs (argv[1], FALSE);

  return (okay ? 0 : 1) ;
}

scorePairs.inc

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <malloc.h>

#include "ArgType.h"

#define LineLen 10000


#ifdef MINIMIZE_WHAT
#define FIND_MIN 1
#else
#define FIND_MIN 0
#endif

// Possibilities for MINIMIZE_WHAT:
// #define MINIMIZE_WHAT pScores->FPFN
// #define MINIMIZE_WHAT pScores->MCCzero



typedef double ValueT;

typedef struct {
  int Neighbour;
  ValueT Score;
  int GroupNum;
} ObjectT;

typedef struct {
  int nRecs;
  int nValsPerRec;
  char ** Names;
  ValueT * Values;
  ObjectT * Objects;
  ValueT minVal;
  ValueT maxVal;
  int nErrNrSrc;
  int nErrNMLo;
  ValueT hiMat;
  ValueT loNonMat;
  ValueT highest;
  int actualPos;
  int actualNeg;
  int truePos;
  int falsePos;
  int trueNeg;
  int falseNeg;
  //ValueT sumFalsePc;
  ValueT FPFN;
  ValueT MCC;
  ValueT MCCzero;
} ScoresT;

/*
  Note that the score between two objects is symmetrical,
  but the "nearest neghbour" relation may not be.
*/


static BOOL GetField (char ** pLine, char * sField, int * fldLen)
// Returns whether another field was found.
// Field may be empty.
{
  //printf ("pLine={%s}", *pLine);
  strcpy (sField, *pLine);
  //printf ("pLine={%s}", *pLine);
  //printf (" sField={%s}", sField);

  if (!*sField) return FALSE;
  char * p = strchr (sField, ',');
  if (p) *p = '\0';
  int nLen = strlen (sField);
  //printf ("nLen=%i ", nLen);
  *pLine += nLen;
  if (p) *pLine += 1;

  *fldLen = nLen;

  return TRUE;
}

static BOOL AllocScores (ScoresT * pScores, int nRecs, int nValues)
{
  printf ("AllocScores\n");

  pScores->nRecs = nRecs;
  pScores->nValsPerRec = nValues;

  pScores->Names = (char **) malloc (nRecs * sizeof (char *));
  if (!pScores->Names) {
    printf ("oom Names");
    return FALSE;
  }

  pScores->Values = (ValueT *) malloc (nRecs * nValues * sizeof (ValueT));
  if (!pScores->Values) {
    free (pScores->Names);
    printf ("oom Values");
    return FALSE;
  }

  pScores->Objects = (ObjectT *) malloc (nRecs * sizeof (ObjectT));
  if (!pScores->Objects) {
    free (pScores->Values);
    free (pScores->Names);
    printf ("oom Objects");
    return FALSE;
  }

  int i;
  for (i=0; i < nRecs; i++) {
    pScores->Names[i] = NULL;
    pScores->Objects[i].Neighbour = -1;
    pScores->Objects[i].Score = -1;
    pScores->Objects[i].GroupNum = -1;
  }
  for (i=0; i < nRecs*nValues; i++) pScores->Values[i] = 0.0;

  return TRUE;
}


static void DeAllocScores (ScoresT * pScores)
{
  int i;
  for (i=0; i < pScores->nRecs; i++) {
    if (pScores->Names[i]) free (pScores->Names[i]);
  }
  free (pScores->Names);
  free (pScores->Values);
  free (pScores->Objects);

  pScores->nRecs = 0;
  pScores->nValsPerRec = 0;
}


static inline ValueT * pointValue (ScoresT * pScores, int recNum, int valNum)
{
  if (recNum >= pScores->nRecs || valNum >= pScores->nValsPerRec) {
    exit (1);
  }

  return &(pScores->Values[recNum*pScores->nValsPerRec + valNum]);
}

static void DumpScores (ScoresT * pScores)
{
  printf ("DumpScores nRecs=%i valsPerRec=%i\n", pScores->nRecs, pScores->nValsPerRec);

  int i, j;
  for (i=0; i < pScores->nRecs; i++) {
    printf ("%i: ", i);
    if (pScores->Names[i]) printf ("[%s]", pScores->Names[i]);

    for (j=0; j < pScores->nValsPerRec; j++) {
      printf (", %g", *pointValue (pScores, i, j));
    }

    printf ("\n");
  }
}

static void DumpNeighbours (ScoresT * pScores)
{
  int i;
  for (i=0; i < pScores->nRecs; i++) {
    printf ("%i: %s, ", i, pScores->Names[i]);
    int neighb = pScores->Objects[i].Neighbour;
    if (neighb >= 0) {
      printf ("%s, %g", pScores->Names[neighb], pScores->Objects[i].Score);
    } else {
      printf ("unknown, unknown");
    }
    printf (",%i", pScores->Objects[i].GroupNum);
    printf ("\n");
  }
}

static BOOL ReadFile (FILE * fin, int * nRecs, int * nValues, ScoresT * pScores)
{
  printf ("ReadFile %s\n", pScores ? "saving" : "not saving");
  BOOL okay = TRUE;
  char NextLine [LineLen];
  char Field [LineLen];

  fseek (fin, 0, SEEK_SET);

  int nLines = 0;
  int expectFlds = -1;

  if (pScores) {
    pScores->minVal = 9e9;
    pScores->maxVal = -9e9;
  }

  while (fgets (NextLine, LineLen, fin)) {
    if (!okay) continue;

    //printf ("%s\n", NextLine);

    char *p = NextLine;

    int nFlds=0;
    int fldLen;
    while (GetField (&p, Field, &fldLen)) {
      // printf ("  [%s]\n", Field);
      if (pScores) {
        if (nFlds == 0) {
          pScores->Names[nLines] = (char *) malloc ((fldLen+1) * sizeof (char));
          if (!pScores->Names[nLines]) {
            okay = FALSE;
          } else {
            strcpy (pScores->Names[nLines], Field);
          }

        } else {
          ValueT * pv = pointValue (pScores, nLines, nFlds-1);
          *pv = atof (Field);
//printf ("*pv=%g ", *pv);

          if (pScores->minVal > *pv) pScores->minVal = *pv;
          if (pScores->maxVal < *pv) pScores->maxVal = *pv;
        }
      }
      nFlds++;
    }
    if (expectFlds == -1) expectFlds = nFlds;
    else if (nFlds != expectFlds) {
      okay = FALSE;
    }
    nLines++;
  }

  if (expectFlds == -1) okay = FALSE;

  *nRecs = nLines;
  *nValues = expectFlds - 1;

  return okay;
}

/*
  Infile is a CSV text file without quotes, no header, one object (eg image) per line.
  First field is a name (typically an image filename).
  This is followed by (n) floating-point values.
  All objects must have same number of values.

  Field separators are commas, with no spaces.


  Works by calculating score for every possible pair of objects.
    score = RMS(difference)
  where
    difference is a value in one object minus the corresponding value in the other object.
    RMS is root-mean-square
*/

static inline ValueT CompareRecs (ScoresT * pScores, int r0, int r1)
{
  ValueT v = 0;

  ValueT * pv0 =  &(pScores->Values[r0*pScores->nValsPerRec]);
  ValueT * pv1 =  &(pScores->Values[r1*pScores->nValsPerRec]);

  int i;
  for (i=0; i < pScores->nValsPerRec; i++) {

#if PROP_DIFF==0
    ValueT diff = (*pv1 - *pv0);
#endif

#if PROP_DIFF==1

// Optimum MIN_DIV was found by trial and error.
#define MIN_DIV 0.25

    ValueT div = *pv1 + *pv0;
    if (abs(div) < MIN_DIV) div = (div < 0)? -MIN_DIV : +MIN_DIV;

    ValueT diff = (*pv1 - *pv0) / div;

/*
    ValueT v0 = *pv0;
    ValueT v1 = *pv1;

    if (v0 < 0) v0 = 0;
    if (v1 < 0) v1 = 0;
    div = v1 + v0;
    diff = (v1 - v0);
*/

#endif

    v += diff*diff;
    pv0++;
    pv1++;
  }
  return sqrt (v / pScores->nValsPerRec);
}

static inline void FindBestScore (ScoresT * pScores, int r0)
{
  printf ("FindBestScore\n");

  ValueT BestScore = 9e+9;
  int BestRec = -1;

  int i;
  for (i=0; i < pScores->nRecs; i++) {
    if (i != r0) {
      ValueT vt = CompareRecs (pScores, r0, i);
      if (BestScore > vt) {
        BestScore = vt;
        BestRec = i;
      }
    }
  }
  if (BestRec == -1) {
    exit (1);
  }
  printf ("%s,%s,%g\n", pScores->Names[r0], pScores->Names[BestRec], BestScore);

  pScores->Objects[r0].Neighbour = BestRec;
  pScores->Objects[r0].Score = BestScore;
}


static inline BOOL IsSourceName (ScoresT * pScores, int r1)
{
  char *p = strchr (pScores->Names[r1], '_');

  if (p) return FALSE;
  else   return TRUE;
}

static inline BOOL DoNamesMatch (ScoresT * pScores, int r0, int r1)
{
#define MAX_LEN 100

  char name0[MAX_LEN];
  char name1[MAX_LEN];

  strncpy (name0, pScores->Names[r0], MAX_LEN);
  strncpy (name1, pScores->Names[r1], MAX_LEN);

  char * p;
  p = strchr (name0, '.');
  if (p) *p = '\0';
  p = strchr (name1, '.');
  if (p) *p = '\0';
  p = strchr (name0, '_');
  if (p) *p = '\0';
  p = strchr (name1, '_');
  if (p) *p = '\0';
  p = strchr (name0, '-');
  if (p) *p = '\0';
  p = strchr (name1, '-');
  if (p) *p = '\0';

  return (strcmp (name0, name1) == 0);
}

static inline void FindFourScores (ScoresT * pScores, int r0)
// For record r0, finds four scores:
// highest and lowest with matching name, and highest and lowest with non-matching name.
// "Name matching" is up to first underscore.
{
  printf ("FindFourScores\n");

  ValueT MatchHi = -1, MatchLo = 9e+9, NonMatchHi = -1, NonMatchLo = 9e+9;
  int rmlo=-1, rmhi=-1, rnmlo=-1, rnmhi=-1;

  ValueT SrcScore = 9e+9;
  int nSrc = -1;

  int i;
  for (i=0; i < pScores->nRecs; i++) {

    if (i != r0) {
      ValueT vt = CompareRecs (pScores, r0, i);

      if (IsSourceName (pScores, i)) {
        if (SrcScore > vt) {
          SrcScore = vt;
          nSrc = i;
        }
      }

      if (pScores->highest < vt) pScores->highest = vt;

      if (DoNamesMatch (pScores, r0, i)) {
        if (MatchHi < vt) {
          MatchHi = vt;
          rmhi = i;
        }
        if (MatchLo > vt) {
          MatchLo = vt;
          rmlo = i;
        }
        if (pScores->hiMat < vt) pScores->hiMat = vt;
      } else {
        if (NonMatchHi < vt) {
          NonMatchHi = vt;
          rnmhi = i;
        }
        if (NonMatchLo > vt) {
          NonMatchLo = vt;
          rnmlo = i;
        }
        if (pScores->loNonMat > vt) pScores->loNonMat = vt;
      }
    }
  }
  char * errNrSrc = "";
  char * errNMLo = "";

  if (nSrc < 0) {
    printf ("nSrc < 0\n");
    exit (1);
  }


  if (!DoNamesMatch (pScores, r0, nSrc)) {
    errNrSrc = "**";
    pScores->nErrNrSrc++;
  }

  if (NonMatchLo <= MatchLo) {
    errNMLo = "**";
    pScores->nErrNMLo++;
  }


  printf ("%s\n", pScores->Names[r0]);
  if (rmlo  >= 0)  printf ("  MatLo %s %g\n", pScores->Names[rmlo], MatchLo);
  if (rmhi  >= 0)  printf ("  MatHi %s %g\n", pScores->Names[rmhi], MatchHi);
  if (rnmlo >= 0) printf ("  NonLo %s %g %s\n", pScores->Names[rnmlo], NonMatchLo, errNMLo);
  if (rnmhi >= 0) printf ("  NonHi %s %g\n", pScores->Names[rnmhi], NonMatchHi);
  if (nSrc  >= 0) printf ("  NrSrc %s %g %s\n", pScores->Names[nSrc], SrcScore, errNrSrc);
}

static void FindBestScores (ScoresT * pScores)
{
  int i;
  for (i=0; i < pScores->nRecs; i++) {
    //FindBestScore (pScores, i);
    FindFourScores (pScores, i);
  }
}

static void FindGroup (ScoresT * pScores, int r0)
{
  // Visit all the objects to find the lowest object number.

  //printf ("\nFindGroup %i: ", r0);

  int lowest = 99999;
  int n = r0;
  while (n != lowest && pScores->Objects[n].GroupNum != lowest) {
    if (lowest > n) lowest = n;
    if (pScores->Objects[n].GroupNum < 0) {
      pScores->Objects[n].GroupNum = lowest;
    } else if (pScores->Objects[n].GroupNum > lowest) {
      pScores->Objects[n].GroupNum = lowest;
    } else if (pScores->Objects[n].GroupNum < lowest) {
      lowest = pScores->Objects[n].GroupNum;
    }
    n = pScores->Objects[n].Neighbour;
    //printf ("low=%i n=%i  ", lowest, n);
  }

  n = r0;
  while (pScores->Objects[n].GroupNum != lowest) {
    pScores->Objects[n].GroupNum = lowest;
    n = pScores->Objects[n].Neighbour;
  }
}

static void FindGroups (ScoresT * pScores)
{
  int i;
  for (i=0; i < pScores->nRecs; i++) {
    if (pScores->Objects[i].GroupNum < 0) FindGroup (pScores, i);
  }
}

static void FindMatches (ScoresT * pScores, ValueT threshold)
// Compare every images with every other image.
// If the RMS < threshold, declare them a match.
// From the names, we know if they should actually be a match.
// So we can count false positives and false negatives.
{
  pScores->actualPos = 0;
  pScores->actualNeg = 0;
  pScores->truePos = 0;
  pScores->falsePos = 0;
  pScores->trueNeg = 0;
  pScores->falseNeg = 0;
  //pScores->sumFalsePc = 0.0;

  int i, j;
  for (i=0; i < pScores->nRecs; i++) {
    for (j=0; j < pScores->nRecs; j++) {
      if (i != j) {
        BOOL actualPos = DoNamesMatch (pScores, i, j);
        ValueT vt = CompareRecs (pScores, i, j);
        if (vt < threshold) {
          if (actualPos) {
            pScores->actualPos++;
            pScores->truePos++;
          } else {
            pScores->actualNeg++;
            pScores->falsePos++;
          }
        } else {
          if (actualPos) {
            pScores->actualPos++;
            pScores->falseNeg++;
          } else {
            pScores->actualNeg++;
            pScores->trueNeg++;
          }
        }
      }
    }
  }
  printf ("threshold=%g", threshold);
  printf (" truePos=%i", pScores->truePos);
  printf (" trueNeg=%i", pScores->trueNeg);
  printf (" falsePos=%i", pScores->falsePos);
  printf (" falseNeg=%i", pScores->falseNeg);
  printf (" actualPos=%i", pScores->actualPos);
  printf (" actualNeg=%i", pScores->actualNeg);

/*---
  ValueT fmPc = 100.0 * pScores->falsePos / (ValueT)pScores->actualPos;
  ValueT fmnPc = 100.0 * pScores->falseNeg / (ValueT)pScores->actualNeg;
  pScores->sumFalsePc = fmPc + fmnPc;

  printf ("falsePosPc=%g ", fmPc);
  printf ("falseNegPc=%g ", fmnPc);
  printf ("sumFalsePc=%g ", pScores->sumFalsePc);
---*/

  ValueT FPrate = pScores->falsePos / (ValueT)pScores->actualNeg;
  ValueT FNrate = pScores->falseNeg / (ValueT)pScores->actualPos;
  pScores->FPFN = FPrate + FNrate;

  ValueT d = (ValueT)(pScores->truePos + pScores->falsePos)
           * (ValueT)(pScores->truePos + pScores->falseNeg)
           * (ValueT)(pScores->trueNeg + pScores->falsePos)
           * (ValueT)(pScores->trueNeg + pScores->falseNeg);

  if (d == 0) d = 1;

  ValueT n = ((ValueT)pScores->truePos * (ValueT)pScores->trueNeg)
           - ((ValueT)pScores->falsePos * (ValueT)pScores->falseNeg);

  pScores->MCC = n / sqrt (d);
  pScores->MCCzero = 1 - pScores->MCC;

  printf (" FPrate=%g", FPrate);
  printf (" FNrate=%g", FNrate);

  printf (" FPFN=%g", pScores->FPFN);

  printf (" MCC=%g", pScores->MCC);
  printf (" MCCzero=%g", pScores->MCCzero);

  printf ("\n");
}

static void FindRangeFalseSumPc (ScoresT * pScores, ValueT t0, ValueT t1, int nSteps)
{
  printf ("FindRangeFalseSumPc\n");

  if (t0 > t1) {
    ValueT T = t0;
    t0 = t1;
    t1 = T;
  }
  ValueT tStep = (t1 - t0) / (ValueT)nSteps;

  ValueT t;

  for (t = t0; t <= t1; t += tStep) {
    FindMatches (pScores, t);
  }
}


#if FIND_MIN==1

static inline int sign0 (ValueT v)
{
#define EPSILON 1e-6

  if (fabs(v) < EPSILON) return 0;
  return (v > 0) ? +1 : -1;
}

static void FindThreshForMin (ScoresT * pScores)
{
  printf ("FindThreshForMin\n");

#define epsilon 0.00001

  ValueT t0 = pScores->loNonMat;
  ValueT t3 = pScores->hiMat;

  ValueT t1 = t0 + (t3 - t0) / 3.0;
  ValueT t2 = t0 + (t3 - t0) * 2/3.0;

  FindMatches (pScores, t0);
  ValueT sumPc0 = MINIMIZE_WHAT;

  FindMatches (pScores, t1);
  ValueT sumPc1 = MINIMIZE_WHAT;

  FindMatches (pScores, t2);
  ValueT sumPc2 = MINIMIZE_WHAT;

  FindMatches (pScores, t3);
  ValueT sumPc3 = MINIMIZE_WHAT;

  BOOL done = FALSE;
  int nIter = 0;

  while (!done) {
    printf ("t=%g, %g, %g, %g ", t0, t1, t2, t3);
    printf ("sumPc=%g, %g, %g, %g\n", sumPc0, sumPc1, sumPc2, sumPc3);

/*---
    if (sumPc1 > sumPc0 && sumPc1 > sumPc2) {
      printf ("Bust\n");
      done = TRUE;
    } else if (fabs(t0 - t2) < epsilon) {
      printf ("Done\n");
      done = TRUE;
    } else if (sumPc0 > sumPc2) {
      // Required t is between t1 and t2
      sumPc0 = sumPc1;
      t0 = t1;
    } else if (sumPc0 < sumPc2) {
      // Required t is between t0 and t1
      sumPc2 = sumPc1;
      t2 = t1;
    } else {
      printf ("Huh?\n");
      done = TRUE;
    }
    t1 = (t0 + t2) / 2.0;
---*/

    int s01 = sign0 (sumPc1 - sumPc0);
    int s12 = sign0 (sumPc2 - sumPc1);
    int s23 = sign0 (sumPc3 - sumPc2);

    //printf ("sign %i, %i, %i\n", s01, s12, s23);

    BOOL dropFirst = FALSE;
    BOOL dropLast = FALSE;

    if (s01==0 && s12==0 && s23==0) {
      printf ("Done\n");
      done = TRUE;
    } else if (s01==-1 && s12==-1) {
      dropFirst = TRUE;
    } else if (s23==+1 && s12==+1) {
      dropLast = TRUE;
    }

    if (!done) {
      if (!dropFirst && !dropLast) {
        if (sumPc0 >= sumPc3) dropFirst = TRUE;
        else dropLast = TRUE;
      }

      if (dropFirst) {
        //printf ("dropFirst\n");
        t0 = t1;
        sumPc0 = sumPc1;
      } else if (dropLast) {
        //printf ("dropLast\n");
        t3 = t2;
        sumPc3 = sumPc2;
      }

      t1 = t0 + (t3 - t0) / 3.0;
      t2 = t0 + (t3 - t0) * 2/3.0;

      FindMatches (pScores, t1);
      sumPc1 = MINIMIZE_WHAT;

      FindMatches (pScores, t2);
      sumPc2 = MINIMIZE_WHAT;
    }
    nIter++;
  }
  printf ("nIter=%i\n", nIter);
  printf ("best at threshold=%g MinValue=%g\n", t1, sumPc1);
  FindMatches (pScores, t1);
}
#endif // FIND_MIN==1


static BOOL ScorePairs (char * InFile, BOOL outAll)
// Returns whether okay.
{
  if (!*InFile) {
    printf ("ScorePairs needs InFile");
    return FALSE;
  }

  FILE * fin;

  BOOL okay = TRUE;

  if (strcmp (InFile, "-")==0) {
    fin = stdin;
  } else {
    fin = fopen (InFile, "rt");
  }

  if (fin == NULL) {
    printf ("Can't open input file\n");
    okay = FALSE;
    return okay;
  }

  int nRecs, nValues;
  if (!ReadFile (fin, &nRecs, &nValues, NULL)) {
    okay = FALSE;
    return okay;
  }

  printf ("nRecs=%i  nValues=%i\n", nRecs, nValues);

  ScoresT Scores;
  okay = AllocScores (&Scores, nRecs, nValues);
  if (!okay) return FALSE;

  Scores.minVal = 0;
  Scores.maxVal = 0;

  if (!ReadFile (fin, &nRecs, &nValues, &Scores)) okay = FALSE;

/*
  DumpScores (&Scores);
*/

  Scores.nErrNrSrc = 0;
  Scores.nErrNMLo = 0;
  Scores.hiMat = 0;
  Scores.loNonMat = 99e99;
  Scores.highest = 0;

  FindBestScores (&Scores);
  printf ("nErrNMLo = %i", Scores.nErrNMLo);
  printf (",  nErrNrSrc = %i", Scores.nErrNrSrc);
  printf (",  highest = %g", Scores.highest);
  printf (",  hiMat = %g", Scores.hiMat);
  printf (",  loNonMat = %g\n", Scores.loNonMat);

  printf ("minVal = %g", Scores.minVal);
  printf (", maxVal = %g\n", Scores.maxVal);

/*
  DumpNeighbours (&Scores);

  FindGroups (&Scores);

  DumpNeighbours (&Scores);
*/

  //FindMatches (&Scores, Scores.hiMat);
  //FindMatches (&Scores, (Scores.hiMat+Scores.loNonMat)/2.0);
  //FindMatches (&Scores, Scores.loNonMat);

#if FIND_MIN==0
  FindRangeFalseSumPc (&Scores, Scores.loNonMat, Scores.hiMat, 10);
#endif

#if FIND_MIN==1
  FindThreshForMin (&Scores);
#endif

  DeAllocScores (&Scores);


  fclose (fin);

  return okay;
}


//int main (int argc, char *argv [])
//{
//  BOOL okay = ScorePairs (argv[1], "out.csv", FALSE);
//
//  return (okay ? 0 : 1) ;
//}

All images on this page were created by the commands shown, using:

%IM%identify -version
Version: ImageMagick 6.9.5-3 Q16 x86 2016-07-22 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2015 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Visual C++: 180040629
Features: Cipher DPC Modules OpenMP 
Delegates (built-in): bzlib cairo flif freetype jng jp2 jpeg lcms lqr openexr pangocairo png ps rsvg tiff webp xml zlib

To improve internet download speeds, some images may have been automatically converted (by ImageMagick, of course) from PNG to JPG.

Source file for this web page is phashtest.h1. To re-create this web page, run "procH1 phashtest".


This page, including the images except where shown otherwise, is my copyright. Anyone is permitted to use or adapt any of the code, scripts or images for any purpose, including commercial use.

Anyone is permitted to re-publish this page, but only for non-commercial use.

Anyone is permitted to link to this page, including for commercial use.


Page version v1.0 22-August-2016.

Page created 30-Sep-2016 17:37:08.

Copyright © 2016 Alan Gibson.