snibgo's ImageMagick pages

Perceptual hash tests

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

• References
• Phash calculation
• Objectives
• Variables
• Methodology
• Proportional difference
• Ranges of values
• Conclusions
• Recommendations for ImageMagick internals
• Future
• Scripts
• Programs

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

References

• Perceptual Hashing for Color Images Using Invariant Moments, Zhenjun Tang, Yumin Dai and Xianquan Zhang, 2011.
• Tests Of Perceptual Hash (PHASH) Compare Metric, Fred Weinhaus.

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% 

@pht_dataset1dir.lis,h

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 

@pht_phtest.lis,h

Colorspaces

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

echo %COLSP3% 

@pht_colsp3.lis

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% 

@pht_colsp_no_0.lis

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

echo %COLSPACES% 

@pht_colspaces.lis

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

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

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

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

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

In the lists above, what colorspaces occur most frequently?

call %PICTBAT%phUnion pht_popcolsp.csv
call csv2tabNH pht_popcolsp

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.

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.

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.
@rem
@rem Updated:
@rem   3-September-2022 for IM v7.
@rem

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 %IMG7%magick ^
    %%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 %IMG7%magick ^
    %%F ^
    -gamma %%V ^
    -set gamma 0.454545 ^
    +depth ^
    %OUTDIR%\!BASE!_gam_%%V.%EXT%

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

  for %%V in (5 10 25 50 75 100) do %IMG7%magick ^
    %%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 %IMG7%magick ^
    -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 %IMG7%magick ^
    -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 %IMG7%magick ^
    %%F ^
    -resize %%V%% ^
    %OUTDIR%\!BASE!_res_%%V.%EXT%

  for %%V in (5 15 30 45 90) do %IMG7%magick ^
    %%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 %IMG7%magick ^
    %%F ^
    -modulate 100,100,%%V ^
    %OUTDIR%\!BASE!_hue_%%V.%EXT%

  for %%V in (70 80 90) do %IMG7%magick ^
    %%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%

  %IMG7%magick ^
    %%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.
@rem
@rem Updated:
@rem   3-September-2022 for IM v7.
@rem

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 (`%IMG7%magick ^
      %%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) ;
//}

%IMG7%magick -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

Page version v1.0 22-August-2016.

Page created 03-Sep-2022 05:41:19.

Copyright © 2022 Alan Gibson.