DM Script, LogPolar image registration

//////////////////////////////////////////////////////////////////////////
// SCRIPT FOR DIGITAL MICROGRAPH, Gatan Inc.
//////////////////////////////////////////////////////////////////////////
// Register rotated, scaled, subimage by LogPolar & PhaseCorrelation    //
//////////////////////////////////////////////////////////////////////////
// Description:                                                         //
//  This class implements the suggested alogrithm from the publication: //
//                                                                      //
//  "Image Registration Using Log Polar Transform and                   //
//   Phase Correlation to Recover Higher Scale"                         //
//   by Jignesh N Sarvaiya; Suprava Patnaik; Kajal Kothari;             //
//   in JOURNAL OF PATTERN RECOGNITION RESEARCH 7 (2012) 90-105         //
//                                                                      //
//////////////////////////////////////////////////////////////////////////
// Version History:                                                     //
//      2013/06/27: Original first draft.                               //
//////////////////////////////////////////////////////////////////////////
number kDebug = 0

// Constants
number TRUE  = 1
number FALSE = 0

number LOWPASS  = 0
number HIGHPASS = 1

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Registration Class                                                                                                                               //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Class RegisterSubImage : object
{
    ////////////////////////////
    // General Initialization //
    ////////////////////////////

    //RegisterSubImage( object self ) { Result( "Created RegisterSubImage object (ID:"+self.ScriptObjectGetID()+")\n"); }
    //~RegisterSubImage( object self ) { Result( "Deleted RegisterSubImage object (ID:"+self.ScriptObjectGetID()+")\n"); }

    ///////////////////
    // Image Filters //
    ///////////////////
    image Create_ModifiedHanningWindow(object self, number sizeX, number sizeY, number cutoff, number width, number isHighPass )
    {
        image   filter
        string  addon = isHighPass ? "high-pass" : "low-pass"
        number  maxXY = max(sizeX,sizeY)

        filter := RealImage("modHanning Window "+addon+" Filter(cutoff="+cutoff+";width="+width+")",8,sizeX,sizeY)
        cutoff = min(1, max(0,cutoff) )

        if (1 < sizeY)  // Filter for Two-Dimensional case (use iradius)
        {
            image cutoffmask := RealImage("",8,maxXY,maxXY)
            cutoffmask = ( iradius/maxXY < cutoff+width ) ? ( (iradius/maxXY > cutoff) ? (cutoff+width-iradius/maxXY)/width : 1 ) : 0
            cutoffmask = min(cutoffmask,tert(abs(icol-maxXY/2)/maxXY<=0.5-width,1,(0.5-width-abs(icol-maxXY/2)/maxXY)/width+1))
            cutoffmask = min(cutoffmask,tert(abs(irow-maxXY/2)/maxXY<=0.5-width,1,(0.5-width-abs(irow-maxXY/2)/maxXY)/width+1))
            cutoffmask = 1-cos(pi()/2*cutoffmask)**2

            filter = warp( cutoffmask, maxXY * icol/sizeX, maxXY * irow/sizeY )
            filter.ImageSetDimensionCalibration(1, (sizeY+1)/2, (1/sizeY), "", 1 )
            filter.SetColorMode(4)
        }
        else    // Filter for one-Dimensional case (use icol instead of iradius)
        {
            image cutoffmask := RealImage("",8,maxXY,1)
            cutoffmask = abs((icol-sizeX/2)/sizeX)
            cutoffmask = tert(cutoffmask<min(cutoff,0.5-width),1,tert(cutoffmask<min(cutoff+width,0.5),(min(cutoff+width,0.5)-cutoffmask)/width,0))
            cutoffmask = 1-cos(pi()/2*cutoffmask)**2
            filter = cutoffmask
        }
        filter.ImageSetDimensionCalibration(0,(sizeX+1)/2, (1/sizeX), "", 1 )

        if (TRUE == isHighPass)
            filter = 1-filter

        return filter
    }

    ///////////////////////
    // Transformations   //
    ///////////////////////
    image CreateLogPolarTransform( object self, image img, number origin_x, number origin_y, number phi_min, number phi_max, number phi_steps, number rad_max, number rad_steps )
    {
        // Create the LogPolar transformation of the image.
        // i.e. map the specified log-polar coordinates onto the cartesian grid
        //      and interpolate the source image.

        image r := RealImage( "log(r) axis", 8, rad_steps )
        r = ( icol/iwidth  ) * log( rad_max )

        image p := RealImage( "phi axis", 8, phi_steps )
        p = phi_min + ( icol/iwidth   ) * ( phi_max - phi_min )
        image LP := RealImage( "LogPolar", 8, rad_steps, phi_steps )

        LP.ImageSetDimensionCalibration(0, 0, log(rad_max)/rad_steps, "log(pixel)", 0 )
        LP.ImageSetDimensionCalibration(1, phi_min, ( phi_max-phi_min)/phi_steps,  "rad", 0 )
        LP = img.warp( origin_x + exp(r[icol,0]) * cos(p[irow,0]), origin_y + exp(r[icol,0])*sin(p[irow,0]) )
        return LP
    }

    image ScaleRotateImage( object self, image img, number scaleFactor, number rotateAngle )
    {
        // Rotated & Scale. Use rotate on the larger image to keep best resolution.
        if ( scaleFactor > 1 )
        {
            number sx,sy
            img.getsize(sx,sy)
            image scaled := RealImage( "scaled",4,sx*scaleFactor,sy*scaleFactor)
            scaled = img.warp( icol/scaleFactor,irow/scaleFactor )
            image rotated := scaled.rotate( rotateAngle )
            return rotated
        }
        else
        {
            image rotated := img.rotate( rotateAngle )
            number sx,sy
            rotated.getsize(sx,sy)
            image scaled := RealImage( "scaled",4,sx*scaleFactor,sy*scaleFactor)
            scaled = rotated.warp( icol/scaleFactor,irow/scaleFactor )
            return scaled;
        }
    }

    image ZeroPadImage( object self, image ref, image img )
    {
        number sx,sy,rx,ry
        ref.GetSize(rx,ry)
        img.GetSize(sx,sy)
        image out := RealImage( "zero-padded", 4, max( sx, rx ), max( sy, ry ) )
        out[icol,irow] = img
        return out
    }

    image ZeroPadImage( object self, image ref, image img, number shftX, number shftY )
    {
        number sx,sy,rx,ry
        ref.GetSize(rx,ry)
        img.GetSize(sx,sy)
        image out := RealImage( "zero-padded", 4, rx, ry )
        out = img[icol-shftX,irow-shftY] * ( ( icol<shftX) || (irow<shftY) || (icol>sx+shftX) || (irow>sy+shftY)? 0 : 1 )
        return out
    }


    //////////////////
    // Correlations //
    //////////////////
    void FindPhaseCorrelationMaximum( object self, image ref, image img, number &maxX, number &maxY, number &maxValue, number applyFilterOption )
    {
        // Both ref & img have to be 2D real-images of same size.
        number sx,sy
        ref.GetSize(sx,sy)
        complexImage power
        if ( TRUE == applyFilterOption )
        {
            image filter := self.Create_ModifiedHanningWindow(sx,sy,0.5,0.2,0)
            power = FFT(real(ref)) * conjugate( FFT(real(img)*filter) )
        }
        else
            power = FFT(real(ref)) * conjugate( FFT(real(img)) )

        power *= 1/modulus(power)
        maxValue = max( modulus(iFFT(power)), maxX, maxY )
        maxX -= (maxX > sx/2) ? sx : 0
        maxY -= (maxY > sy/2) ? sy : 0
    }

    void FindScaleAndRotation( object self, image ref, image sense, number &scale, number &rotation, number &shiftX, number &shiftY )
    {
        number preScale = scale
        number preRotate = rotation
        image senseImage
        number rx, ry
        ref.GetSize(rx,ry)

        if (( 1 == preScale ) && ( 0 == preRotate ) )
            senseImage := sense
        else
            senseImage := self.ScaleRotateImage( sense, 1/preScale, rotation )

        senseImage := self.ZeroPadImage( ref, senseImage, shiftX, shiftY )

        image frequency_filter := self.Create_ModifiedHanningWindow(rx,ry,0.01,0.1,TRUE)        // preferably use high-frequency info
        image space_filter := self.Create_ModifiedHanningWindow(rx,ry,0.5,0.2,FALSE)            // filter before FFT to remove 'cross' artefacts

        image F_Ref  := modulus( FFT( real(ref*space_filter) ) )*frequency_filter
        image F_Sense  := modulus( FFT( real(senseImage*space_filter) ) )*frequency_filter

        number angular_min = 0;
        number angular_max = 2*Pi();
        number angular_step = 760;
        number rad_max = maximum(rx/2,ry/2 )
        number rad_step= rad_max*2

        image ref_PL   := self.CreateLogPolarTransform( F_Ref  , rx/2, ry/2, angular_min, angular_max, angular_step, rad_max, rad_step )
        image sense_PL := self.CreateLogPolarTransform( F_Sense, rx/2, ry/2, angular_min, angular_max, angular_step, rad_max, rad_step )
        if ( 1 == kDebug )
        {
            ref_PL.showImage()
            sense_PL.showImage()
            exit(0)
        }
        number maxX, maxY, maxV
        self.FindPhaseCorrelationMaximum( ref_PL, sense_PL, maxX, maxY, maxV, TRUE )

        rotation = (angular_max - angular_min) / angular_step * (-maxY)
        scale    = exp( log(rad_step)/rad_step ) ** (maxX)
        scale*=prescale
        rotation+=prerotate

        //result("\n\n\t Found shifts :x="+maxX+" y="+maxY)
        //result("\n\t Found scale: "+scale+" ("+prescale+")")
        //result("\n\t Found rotation : "+rotation/Pi()*180+" degree, or "+(rotation/Pi()*180-180)+" degree ("+preRotate/Pi()*180+")\n\n")
    }

    void FindTranslationScaleAndRotation( object self, image ref, image sense, number &scale, number &rotation, number &ShiftX, number &ShiftY )
    {
        // First determine scale & rotation in a translational-invariate manner
        result("\nInitial scale    : "+scale+"\n")
        result("Initial rotation : "+rotation/Pi()*180+" [degree]\n")
        result("Initial shifts XY: "+ShiftX+"/"+ShiftY+" [pixel]\n\n")

        OpenAndSetProgressWindow( "Find T&S&R","get rought S&R","")
        self.FindScaleAndRotation( ref, sense, scale, rotation, shiftX, shiftY );

        OpenAndSetProgressWindow( "Find T&S&R","find T","S="+Format(scale,"%3.2f")+" R="+Format( rotation*180/Pi(), "%5.1f") )
        image senseImage1 := self.ScaleRotateImage( sense, 1/scale, rotation )
        senseImage1 := self.ZeroPadImage( ref, senseImage1, shiftX, shiftY )

        image senseImage2 := self.ScaleRotateImage( sense, 1/scale, rotation-Pi() )
        senseImage2 := self.ZeroPadImage( ref, senseImage2, shiftX, shiftY)

        number maxX1, maxY1, maxV1
        self.FindPhaseCorrelationMaximum( ref, senseImage1, maxX1, maxY1, maxV1, TRUE )
        number maxX2, maxY2, maxV2
        self.FindPhaseCorrelationMaximum( ref, senseImage2, maxX2, maxY2, maxV2, TRUE )
        ShiftX += (maxV1>maxV2) ? maxX1 : maxX2
        ShiftY += (maxV1>maxV2) ? maxY1 : maxY2
        rotation = (maxV1>maxV2) ? rotation : rotation-Pi()
        rotation += (rotation<0) ? 2*Pi() : 0
        result("Found scale      : "+scale+"\n")
        result("Found rotation   : "+rotation/Pi()*180+" [degree]\n")
        result("Found shifts XY  : "+ShiftX+"/"+ShiftY+" [pixel]\n")
        result("Phase correlation: "+max(maxV1,maxV2)+"\n")
        OpenAndSetProgressWindow( "","","")
    }

    ///////////////////////
    // Results & Display //
    ///////////////////////
    image RGBOverlay( object self, image ref, image sense, number scale, number rotation, number shftX, number shftY, number &getCorMax )
    {
        image sen := self.ScaleRotateImage( sense, 1/scale, rotation )
        image cut := self.ZeroPadImage( ref, sen, shftX, shftY )
        if ( TRUE == getCorMax )
        {
            image ref_cut = ( cut == 0 ) ? 0 : ref
            image crossCorr := CrossCorrelate( cut, ref_cut )
            getCorMax = max(crossCorr)
            Result( "\n Correlation: " + getCorMax  +"\n")
        }

        return RGB(ref/max(ref)*255,cut/max(cut)*255,0)
    }

    number GetCorrelationCoefficient( object self, image ref, image sense,  number scale, number rotation, number shftX, number shftY )
    {
        image sen := self.ScaleRotateImage( sense, 1/scale, rotation )
        image sen_cut := self.ZeroPadImage( ref, sen, shftX, shftY )
        image ref_cut = ( sen_cut == 0 ) ? 0 : ref
        image crossCorr := CrossCorrelate( sen_cut, ref_cut )
        return max(crossCorr)
    }
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// MAIN PROGRAM (Testing RegisterClass)                                                                                                                         //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
{
    kDebug = 0
    image ref := B
    image sen := D
    object ob = Alloc(RegisterSubImage)
    number mx,my,mv
    number scal,rot,shftx,shfty
    // If a rough-aligment has already been made,
    // it can be fed in here to improve results of the algorithm.
    scal    = 4//1.7 //1.7
    rot     = 0//-45 *PI()/180
    shftx   = 0//201
    shfty   = 0//270
    //ob.ScaleRotateImage( sen, 1/scal, rot ).showimage()
    //exit(0)
    number CorMax = TRUE
    ClearResults()
    ob.FindTranslationScaleAndRotation(ref,sen,scal,rot,shftx,shfty )
    ob.RGBOverlay(ref,sen,scal,rot,shftx,shfty,CorMax).ShowImage()
}