snibgo's ImageMagick pages

Perceptual hash thresholds

We experiment with RMS thresholds that classify image pairs as "same" or "different".

• Methodology
• Conclusions
• Recommendations for ImageMagick internals
• Scripts
• Programs

This page continues testing from Perceptual hash tests. (For references, see that page.) That page took a small number of standard images, created many variations ("attacks") of each, and calculated the phash of every image, in combinations of up to four colorspaces. For every colorspace combination, it tested RMSE phash diferences between every pair of images. Two measures of "goodness" were used, according to how many images were closest to an image that wasn't derived from the same standard image, or how many were closest to a standard image that wasn't the actual source for the variation.

On this page, we use measures of false positive and false negative matches at different RMS thresholds to find the "best" threshold for each colorspace combination. The measures are: sum of false rates, and Matthews correlation coefficient. These measures also indicate the "goodness" of each colorspace combination, so we get a ranking order of the combinations.

The pages have different criteria for goodness, so generate different ranking orders.

Methodology

The page Perceptual hash tests has calculated phash values for every image, in a large number of colorspace combinations. From these numbers, we can calculate RMS differences between any pair of images; this number is zero when the images are identical, and larger numbers mean the images are more different. At each combination, it has found the lowest RMS difference between images that were variations of different standard images (loNonMat), and the largest RMS difference between images that were variations of the same standard images (hiMat).

For this image set, loNonMat < hiMat.

If we choose a threshold less than loNonMat, all the images from different sources will score above the threshold so will correctly be considered different; there will be no false positives. However, there will be many false negatives.

If we choose a threshold greater than hiMat, all the images from the same source will score below the threshold so will correctly be considered the same; there will be no false negatives. However, there will be many false positives.

If we choose a threshold between those limits, we get false positives and false negatives.

The script setPhEnv.bat sets up the environment:

call %PICTBAT%setPhEnv

This page needs Inkscape:

call %PICTBAT%setInkPath

We want the "optimum" threshold for the RMS phash. But what defines optimum? What do we want to optimise?

True/false positive/negative

References:

• Wikipedia: Detection error tradeoff.
• Wikipedia: Receiver operating characteristic.
• Wikipedia: Confusion matrix.
• Wikipedia: Matthews correlation coefficient.
• Evaluation: From Precision, Recall and F-Factor to ROC, Informedness, Markedness & Correlation, David M W Powers, 2007.

Of all the possible pairs, we know which ones should match (eg every Lena variation should match every other Lena variation), and which ones shouldn't (eg a Lena variation shouldn't match a Barbara variation).

If the RMS difference betwen the pair is below the given threshold, we call it a positive match, and we also know if this is a true positive or false positive. Similarly when the RMS difference is above the threshold, we call it a negative match, we know if this is a true negative or false negative.

With this dataset, no threshold gives zero false positives and zero false negatives. We can only tradeoff one against the other.

We will draw Detection Error Tradeoff (DET) curves. This will map the false negative rate against the false positive rate, for various values of RMS threshold.

gensvg /itrueFalsePosNeg.s1 /P

If we have R records, the number of pairs is nPairs = R*(R-1)/2. In the diagram, the total rectangle represents nPairs.

Of these nPairs pairs, actualPos pairs will match and actualNeg = (nPairs - actualPos) will not match. In the diagram, the green shape represents actualPos, and everything outside the green shape represents actualNeg.

Each test run will consider p pairs to match and q = (nPairs - p) pairs to not match. In the diagram, the red shape represents p.

If the method was perfect for the dataset, the red shape would coincide with the green shape, so falseNeg and falsePos would both be zero.

Of the p pairs that are considered to match, truePos are true matches and falsePos = (p - truePos) are false matches.

Of the q pairs that are considered to not match, trueNeg are true non-matches and falseNeg = (q - trueNeg) are false non-matches.

Hence:

truePos + falseNeg = actualPos

trueNeg + falsePos = actualNeg

Instead of using integer counts of pairs, we divide to make false positive and false negative rates that range from 0.0 (a perfect score) to 1.0 (as bad as it could be.)

FPrate = falsePos / actualNeg = falsePos / (trueNeg + falsePos)

FNrate = falseNeg / actualPos = falseNeg / (truePos + falseNeg)

We will plot the FNrate (y-axis) against FPrate (x-axis) for a number of thresholds.

call %PICTBAT%phpairs 21_HSB
call %PICTBAT%phpairs 21_YDbDr
call %PICTBAT%phpairs 42_HSB_YDbDr
call %PICTBAT%phpairs 42_HSB_xyY

At any FPrate, the YDbDr colorspace gives a better (nearer zero) FNrate than the HSB colorspace. And the combination of two colorspaces HSB+YDrBr is better still. Where many points are plotted close together, this is where scorePairsNPD has searched for a minimum of (FPrate + FNrate).

We can plot all the single-colorspaces on one graph, zooming in to a small range.

Okay, but what is the best threshold value? I'm not a statistician. A naive "best" threshold is that which minimizes FPFN = (FPrate + FNrate). Howevever, perhaps the threshold that maximizes the Matthews Correlation Coefficient is more appropriate.

MCC values range from +1 through 0 to -1. MCC = +1 means the method was entirely correct. MCC = 0 means the method was no better than random guesses. MCC = -1 means the method gave entirely the wrong answers. As I prefer to minimize rather than maximize, I use MCCzero = (1 - MCC).

MCC is calculated by:

n = truePos * trueNeg - falsePos * falseNeg

d = (truePos + falsePos) * (truePos + falseNeg) * (trueNeg + falsePos) * (trueNeg + falseNeg)

if d==0 then d = 1

MCC = n / sqrt (d)

This provides alternate measures for the best colorspace combination: it is the colorspace that yields the lowest FPFN or MCCzero.

The program scorePairsMaxMcc.c iterates to find the threshold, for one colorspace combination, that minimizes MCCzero. This takes around six seconds (single-threaded). (The program would be faster if we stored the RMS differences of the 6,859,841,200 comparisons.)

The script phpPairsBestAll.bat repeats this for all 31,930 combinations of colorspace, which takes around 53 hours.

Results

We list the results in decreasing order of MCC, so the colorspace combinations at the top are the best.

call %PICTBAT%phpPairsBestAll phpbest_XX.csv
if ERRORLEVEL 1 goto error

cHead /iphpbest_2.csv /ophpbest_2_15.csv /h16
cHead /iphpbest_3.csv /ophpbest_3_15.csv /h16
cHead /iphpbest_4.csv /ophpbest_4_15.csv /h16

call csv2tab phpbest_1
call csv2tab phpbest_2_15
call csv2tab phpbest_3_15
call csv2tab phpbest_4_15

The results for all of the single colorspaces:

The combinations of two colorspaces, showing just the top 15 out of 435:

The combinations of three colorspaces, showing just the top 15 out of 4060:

The combinations of four colorspaces, showing just the top 15 out of 27,405:

In the lists above, what colorspaces occur most frequently?

call %PICTBAT%phUnionTh phh_popcolsp.csv
call csv2tabNH phh_popcolsp

Conclusions

The method on Perceptual hash tests page measures whether the closest match to a given image is a correct one, and doesn't care about whether more distant matches are also correct.

By contrast, this page measures the goodness of colorspace combination by a calculation from true or false positives or negatives, given a certain RMS threshold. This measure is appropriate when the application is to find images that match and reject those that don't.

Recommendations for ImageMagick internals

My recommendation for IM internals are the same as shown on Perceptual hash tests, except that the recommended colorspaces, when two are used, was given there as HSB+xyY. On ths page, we find that HSB+xyY has an MCC score of 0.783679. The best MCC score, of 0.802692, comes from OHTA+UV. This is not a large difference; when the test is MCC of false positives and negatives, using two colorspaces, the test is not highly sensitive to the choice of colorspaces.

In this test, HSB and xyY appear among the top three colorspaces in the "frequently occur" list of best results. The poll positon is taken by OHTA, which scored notably badly on the Perceptual hash tests page.

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%

phpairs.bat

@rem
@rem Updated:
@rem   15-August-2022 for IM v7.
@rem

rem %1 is colorspace combination eg 63_HSB_IQ_xyY

call %PICTBAT%setPhEnv

set COLSP_COMB=%1

set TITLE=%COLSP_COMB:_=+%
set TITLE=%TITLE:21+=%
set TITLE=%TITLE:42+=%
set TITLE=%TITLE:63+=%
set TITLE=%TITLE:84+=%

set COLSP_CSV=%PH_CSVDIR%\idPH%COLSP_COMB%.csv

set TMP_RESULTS=%TEMP%\php.csv

call %PICTBAT%phpThreshCsv %COLSP_CSV% %TMP_RESULTS%
if ERRORLEVEL 1 exit /B 1

set PLOT=php_%COLSP_COMB%

gnuplot-base ^
  -c %PICTBAT%plotScr.gp ^
  %TMP_RESULTS% %PLOT%.svg 500 500 "2" ^
  "FPrate" "FNrate" ^
  "set key inside top right noautotitle title '%TITLE%';set xrange [0:1];set yrange [0:1];set size square"

rem  "set key inside top right noautotitle title '%TITLE%';set logscale;set xrange [0.001:1];set yrange [0.001:1];set size square"

if ERRORLEVEL 1 exit /B 1

rem set yrange [0:1];

%IMG7%magick ^
  -background None %PLOT%.svg ^
  -trim +repage ^
  -bordercolor None -border 5 ^
  -background #eee -layers flatten ^
  %PLOT%.png

phpThreshCsv.bat

rem From %1, a CSV of PH values for images,
rem makes %2 a CSV with varying thresholds of RMS.

setlocal

set COLSP_CSV=%1

set OUTFILE=%2

set TMP_RESULTS=%TEMP%\ptc.csv

if not exist %COLSP_CSV% (
  echo Can't find %COLSP_CSV%
  exit /B 1
)

scorePairsNPD %COLSP_CSV% >%TMP_RESULTS%

cGrep    /p0 /i%TMP_RESULTS% /o%TMP_RESULTS%  /sFPFN

chStrs   /p0 /i%TMP_RESULTS% /o%TMP_RESULTS% /f", " /t","
chStrs   /p0 /i%TMP_RESULTS% /o%TMP_RESULTS% /f" " /t","

call deEqCols %TMP_RESULTS% %OUTFILE%
if ERRORLEVEL 1 exit /B 1

cProject /p0 /i%OUTFILE% /o%OUTFILE% /h /kFPrate,FNrate,FPFN,MCC,threshold

if ERRORLEVEL 1 exit /B 1

cSort    /p0 /i%OUTFILE% /o%OUTFILE% /h /kFPrate

type %TMP_RESULTS%

endlocal

phpThreshCsvBest.bat

rem From %1, a CSV of PH values for images,
rem makes %2 a CSV with best thresholds of RMS,
rem with colorspace column value %3.

setlocal

set COLSP_CSV=%1

set OUTFILE=%2

set TMP_RESULTS=%TEMP%\ptc.csv

if not exist %COLSP_CSV% (
  echo Can't find %COLSP_CSV%
  exit /B 1
)

scorePairsMaxMcc %COLSP_CSV% >%TMP_RESULTS%

cGrep    /p0 /i%TMP_RESULTS% /o%TMP_RESULTS%  /sbest /t\0

cGrep    /p0 /i%TMP_RESULTS% /o%TMP_RESULTS%  /sFPFN

chStrs   /p0 /i%TMP_RESULTS% /o%TMP_RESULTS% /f", " /t","
chStrs   /p0 /i%TMP_RESULTS% /o%TMP_RESULTS% /f" " /t","
cPrefix  /p0 /i%TMP_RESULTS% /o%TMP_RESULTS% /r",colorspace=%3"

call deEqCols %TMP_RESULTS% %OUTFILE%
if ERRORLEVEL 1 exit /B 1

type %OUTFILE%

endlocal

phpPairsBestAll.bat

rem Creates output CSV files with best results.
rem %1 is name format, must contain XX.
rem XX will be replaced with 1, 2, 3 or 4.

setlocal

call %PICTBAT%setPhEnv

set OUT_FMT=%1

set USE_COLS=threshold,truePos,trueNeg,falsePos,falseNeg,FPrate,FNrate,FPFN,MCC,colorspace

rem ------------------------------------------------
rem One colorspace.

goto skip1

set OUTFILE=%OUT_FMT:XX=1%

del phbest21*.csv 2>nul

for %%C in (%COLSPACES%) do (
  call %PICTBAT%phCombin %%C
  if ERRORLEVEL 1 exit /B 1

  set INCSV=%PH_CSVDIR%\idPh!PH_COMBIN!.csv

  set OUTCSV=phbest!PH_COMBIN!.csv

  call %PICTBAT%phpThreshCsvBest !INCSV! !OUTCSV! !PH_PLUS!
  if ERRORLEVEL 1 exit /B 1
)

cat phbest21*.csv >%OUTFILE%

cSort    /p0 /i%OUTFILE% /o%OUTFILE% /u /r
cProject /p0 /i%OUTFILE% /o%OUTFILE% /h /k%USE_COLS%
cSort    /p0 /i%OUTFILE% /o%OUTFILE% /h /kMCC /r

type %OUTFILE%

:skip1


rem ------------------------------------------------
rem Combinations of two colorspaces.

goto skip2

set OUTFILE=%OUT_FMT:XX=2%

del phbest42*.csv 2>nul

for %%C in (%COLSPACES%) do (
  for %%D in (%COLSPACES%) do (
    if /I %%C LSS %%D (

      call %PICTBAT%phCombin %%C %%D
      if ERRORLEVEL 1 exit /B 1

      set INCSV=%PH_CSVDIR%\idPh!PH_COMBIN!.csv

      set OUTCSV=phbest!PH_COMBIN!.csv

      call %PICTBAT%phpThreshCsvBest !INCSV! !OUTCSV! !PH_PLUS!
      if ERRORLEVEL 1 exit /B 1
    )
  )
)

cat phbest42*.csv >%OUTFILE%

cSort /p0 /i%OUTFILE% /o%OUTFILE% /u /r
cSort /p0 /i%OUTFILE% /o%OUTFILE% /h /kMCCzero

cProject /p0 /i%OUTFILE% /o%OUTFILE% /h /k%USE_COLS%

type %OUTFILE%

:skip2


rem ------------------------------------------------
rem Combinations of three colorspaces.

goto skip3

set OUTFILE=%OUT_FMT:XX=3%

del phbest63*.csv 2>nul

for %%C in (%COLSPACES%) do (
  for %%D in (%COLSPACES%) do (
    if /I %%C LSS %%D (
      for %%E in (%COLSPACES%) do (
        if /I %%D LSS %%E (

          call %PICTBAT%phCombin %%C %%D %%E
          if ERRORLEVEL 1 exit /B 1

          set INCSV=%PH_CSVDIR%\idPh!PH_COMBIN!.csv

          set OUTCSV=phbest!PH_COMBIN!.csv

          call %PICTBAT%phpThreshCsvBest !INCSV! !OUTCSV! !PH_PLUS!
          if ERRORLEVEL 1 exit /B 1
        )
      )
    )
  )
)

cat phbest63*.csv >%OUTFILE%

cSort /p0 /i%OUTFILE% /o%OUTFILE% /u /r
cSort /p0 /i%OUTFILE% /o%OUTFILE% /h /kMCCzero

cProject /p0 /i%OUTFILE% /o%OUTFILE% /h /k%USE_COLS%

type %OUTFILE%

:skip3

rem ------------------------------------------------
rem Combinations of four colorspaces.

rem goto skip4

set OUTFILE=%OUT_FMT:XX=4%

del phbest84*.csv 2>nul

for %%C in (%COLSPACES%) do (
  for %%D in (%COLSPACES%) do (
    if /I %%C LSS %%D (
      for %%E in (%COLSPACES%) do (
        if /I %%D LSS %%E (
          for %%F in (%COLSPACES%) do (
            if /I %%E LSS %%F (

              call %PICTBAT%phCombin %%C %%D %%E %%F
              if ERRORLEVEL 1 exit /B 1

              set INCSV=%PH_CSVDIR%\idPh!PH_COMBIN!.csv

              set OUTCSV=phbest!PH_COMBIN!.csv

              call %PICTBAT%phpThreshCsvBest !INCSV! !OUTCSV! !PH_PLUS!
              if ERRORLEVEL 1 exit /B 1
            )
          )
        )
      )
    )
  )
)

cat phbest84*.csv >%OUTFILE%

cSort /p0 /i%OUTFILE% /o%OUTFILE% /u /r
cSort /p0 /i%OUTFILE% /o%OUTFILE% /h /kMCCzero

cProject /p0 /i%OUTFILE% /o%OUTFILE% /h /k%USE_COLS%

type %OUTFILE%

:skip4


endlocal

Programs

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

scorePairsMaxMcc.c

// Whether to use proportional difference.
#define PROP_DIFF 0
#define MINIMIZE_WHAT pScores->MCCzero

#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) ;
//}

%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

Page version v1.0 1-September-2016.

Page created 15-Aug-2022 16:28:32.

Copyright © 2022 Alan Gibson.