Guest User

DM Script, LogPolar image registration

a guest
Aug 1st, 2018
89
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. //////////////////////////////////////////////////////////////////////////
  2. // SCRIPT FOR DIGITAL MICROGRAPH, Gatan Inc.
  3. //////////////////////////////////////////////////////////////////////////
  4. // Register rotated, scaled, subimage by LogPolar & PhaseCorrelation    //
  5. //////////////////////////////////////////////////////////////////////////
  6. // Description:                                                         //
  7. //  This class implements the suggested alogrithm from the publication: //
  8. //                                                                      //                         
  9. //  "Image Registration Using Log Polar Transform and                   //
  10. //   Phase Correlation to Recover Higher Scale"                         //
  11. //   by Jignesh N Sarvaiya; Suprava Patnaik; Kajal Kothari;             //
  12. //   in JOURNAL OF PATTERN RECOGNITION RESEARCH 7 (2012) 90-105         //
  13. //                                                                      //
  14. //////////////////////////////////////////////////////////////////////////
  15. // Version History:                                                     //
  16. //      2013/06/27: Original first draft.                               //
  17. //////////////////////////////////////////////////////////////////////////
  18. number kDebug = 0
  19.  
  20. // Constants
  21. number TRUE  = 1
  22. number FALSE = 0
  23.  
  24. number LOWPASS  = 0
  25. number HIGHPASS = 1
  26.  
  27. //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  28. // Registration Class                                                                                                                               //
  29. //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  30. Class RegisterSubImage : object
  31. {
  32.     ////////////////////////////
  33.     // General Initialization //
  34.     ////////////////////////////
  35.    
  36.     //RegisterSubImage( object self ) { Result( "Created RegisterSubImage object (ID:"+self.ScriptObjectGetID()+")\n"); }
  37.     //~RegisterSubImage( object self ) { Result( "Deleted RegisterSubImage object (ID:"+self.ScriptObjectGetID()+")\n"); }
  38.  
  39.     ///////////////////
  40.     // Image Filters //
  41.     ///////////////////
  42.     image Create_ModifiedHanningWindow(object self, number sizeX, number sizeY, number cutoff, number width, number isHighPass )
  43.     {
  44.         image   filter
  45.         string  addon = isHighPass ? "high-pass" : "low-pass"
  46.         number  maxXY = max(sizeX,sizeY)
  47.            
  48.         filter := RealImage("modHanning Window "+addon+" Filter(cutoff="+cutoff+";width="+width+")",8,sizeX,sizeY)
  49.         cutoff = min(1, max(0,cutoff) )
  50.        
  51.         if (1 < sizeY)  // Filter for Two-Dimensional case (use iradius)
  52.         {
  53.             image cutoffmask := RealImage("",8,maxXY,maxXY)
  54.             cutoffmask = ( iradius/maxXY < cutoff+width ) ? ( (iradius/maxXY > cutoff) ? (cutoff+width-iradius/maxXY)/width : 1 ) : 0
  55.             cutoffmask = min(cutoffmask,tert(abs(icol-maxXY/2)/maxXY<=0.5-width,1,(0.5-width-abs(icol-maxXY/2)/maxXY)/width+1))
  56.             cutoffmask = min(cutoffmask,tert(abs(irow-maxXY/2)/maxXY<=0.5-width,1,(0.5-width-abs(irow-maxXY/2)/maxXY)/width+1))
  57.             cutoffmask = 1-cos(pi()/2*cutoffmask)**2
  58.  
  59.             filter = warp( cutoffmask, maxXY * icol/sizeX, maxXY * irow/sizeY )
  60.             filter.ImageSetDimensionCalibration(1, (sizeY+1)/2, (1/sizeY), "", 1 )
  61.             filter.SetColorMode(4)
  62.         }
  63.         else    // Filter for one-Dimensional case (use icol instead of iradius)
  64.         {
  65.             image cutoffmask := RealImage("",8,maxXY,1)
  66.             cutoffmask = abs((icol-sizeX/2)/sizeX) 
  67.             cutoffmask = tert(cutoffmask<min(cutoff,0.5-width),1,tert(cutoffmask<min(cutoff+width,0.5),(min(cutoff+width,0.5)-cutoffmask)/width,0))
  68.             cutoffmask = 1-cos(pi()/2*cutoffmask)**2
  69.             filter = cutoffmask
  70.         }
  71.         filter.ImageSetDimensionCalibration(0,(sizeX+1)/2, (1/sizeX), "", 1 )
  72.        
  73.         if (TRUE == isHighPass)
  74.             filter = 1-filter
  75.            
  76.         return filter
  77.     }
  78.    
  79.     ///////////////////////
  80.     // Transformations   //
  81.     ///////////////////////
  82.     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 )
  83.     {
  84.         // Create the LogPolar transformation of the image.
  85.         // i.e. map the specified log-polar coordinates onto the cartesian grid
  86.         //      and interpolate the source image.
  87.        
  88.         image r := RealImage( "log(r) axis", 8, rad_steps )
  89.         r = ( icol/iwidth  ) * log( rad_max )  
  90.        
  91.         image p := RealImage( "phi axis", 8, phi_steps )
  92.         p = phi_min + ( icol/iwidth   ) * ( phi_max - phi_min )
  93.         image LP := RealImage( "LogPolar", 8, rad_steps, phi_steps )
  94.        
  95.         LP.ImageSetDimensionCalibration(0, 0, log(rad_max)/rad_steps, "log(pixel)", 0 )
  96.         LP.ImageSetDimensionCalibration(1, phi_min, ( phi_max-phi_min)/phi_steps,  "rad", 0 )
  97.         LP = img.warp( origin_x + exp(r[icol,0]) * cos(p[irow,0]), origin_y + exp(r[icol,0])*sin(p[irow,0]) )
  98.         return LP
  99.     }
  100.    
  101.     image ScaleRotateImage( object self, image img, number scaleFactor, number rotateAngle )
  102.     {
  103.         // Rotated & Scale. Use rotate on the larger image to keep best resolution.
  104.         if ( scaleFactor > 1 )
  105.         {
  106.             number sx,sy
  107.             img.getsize(sx,sy)
  108.             image scaled := RealImage( "scaled",4,sx*scaleFactor,sy*scaleFactor)
  109.             scaled = img.warp( icol/scaleFactor,irow/scaleFactor )
  110.             image rotated := scaled.rotate( rotateAngle )
  111.             return rotated
  112.         }
  113.         else
  114.         {
  115.             image rotated := img.rotate( rotateAngle )
  116.             number sx,sy
  117.             rotated.getsize(sx,sy)
  118.             image scaled := RealImage( "scaled",4,sx*scaleFactor,sy*scaleFactor)
  119.             scaled = rotated.warp( icol/scaleFactor,irow/scaleFactor )
  120.             return scaled;
  121.         }  
  122.     }
  123.    
  124.     image ZeroPadImage( object self, image ref, image img )
  125.     {
  126.         number sx,sy,rx,ry
  127.         ref.GetSize(rx,ry)
  128.         img.GetSize(sx,sy)
  129.         image out := RealImage( "zero-padded", 4, max( sx, rx ), max( sy, ry ) )
  130.         out[icol,irow] = img
  131.         return out
  132.     }
  133.    
  134.     image ZeroPadImage( object self, image ref, image img, number shftX, number shftY )
  135.     {
  136.         number sx,sy,rx,ry
  137.         ref.GetSize(rx,ry)
  138.         img.GetSize(sx,sy)
  139.         image out := RealImage( "zero-padded", 4, rx, ry )
  140.         out = img[icol-shftX,irow-shftY] * ( ( icol<shftX) || (irow<shftY) || (icol>sx+shftX) || (irow>sy+shftY)? 0 : 1 )
  141.         return out
  142.     }
  143.    
  144.    
  145.     //////////////////
  146.     // Correlations //
  147.     //////////////////
  148.     void FindPhaseCorrelationMaximum( object self, image ref, image img, number &maxX, number &maxY, number &maxValue, number applyFilterOption )
  149.     {
  150.         // Both ref & img have to be 2D real-images of same size.
  151.         number sx,sy
  152.         ref.GetSize(sx,sy)
  153.         complexImage power
  154.         if ( TRUE == applyFilterOption )
  155.         {
  156.             image filter := self.Create_ModifiedHanningWindow(sx,sy,0.5,0.2,0)
  157.             power = FFT(real(ref)) * conjugate( FFT(real(img)*filter) )
  158.         }
  159.         else
  160.             power = FFT(real(ref)) * conjugate( FFT(real(img)) )
  161.            
  162.         power *= 1/modulus(power)
  163.         maxValue = max( modulus(iFFT(power)), maxX, maxY )
  164.         maxX -= (maxX > sx/2) ? sx : 0
  165.         maxY -= (maxY > sy/2) ? sy : 0
  166.     }
  167.    
  168.     void FindScaleAndRotation( object self, image ref, image sense, number &scale, number &rotation, number &shiftX, number &shiftY )
  169.     {
  170.         number preScale = scale
  171.         number preRotate = rotation
  172.         image senseImage
  173.         number rx, ry
  174.         ref.GetSize(rx,ry)
  175.    
  176.         if (( 1 == preScale ) && ( 0 == preRotate ) )
  177.             senseImage := sense
  178.         else
  179.             senseImage := self.ScaleRotateImage( sense, 1/preScale, rotation )
  180.  
  181.         senseImage := self.ZeroPadImage( ref, senseImage, shiftX, shiftY )
  182.    
  183.         image frequency_filter := self.Create_ModifiedHanningWindow(rx,ry,0.01,0.1,TRUE)        // preferably use high-frequency info
  184.         image space_filter := self.Create_ModifiedHanningWindow(rx,ry,0.5,0.2,FALSE)            // filter before FFT to remove 'cross' artefacts
  185.            
  186.         image F_Ref  := modulus( FFT( real(ref*space_filter) ) )*frequency_filter
  187.         image F_Sense  := modulus( FFT( real(senseImage*space_filter) ) )*frequency_filter
  188.        
  189.         number angular_min = 0;
  190.         number angular_max = 2*Pi();
  191.         number angular_step = 760;
  192.         number rad_max = maximum(rx/2,ry/2 )
  193.         number rad_step= rad_max*2
  194.        
  195.         image ref_PL   := self.CreateLogPolarTransform( F_Ref  , rx/2, ry/2, angular_min, angular_max, angular_step, rad_max, rad_step )
  196.         image sense_PL := self.CreateLogPolarTransform( F_Sense, rx/2, ry/2, angular_min, angular_max, angular_step, rad_max, rad_step )
  197.         if ( 1 == kDebug )
  198.         {
  199.             ref_PL.showImage()
  200.             sense_PL.showImage()
  201.             exit(0)
  202.         }
  203.         number maxX, maxY, maxV
  204.         self.FindPhaseCorrelationMaximum( ref_PL, sense_PL, maxX, maxY, maxV, TRUE )
  205.        
  206.         rotation = (angular_max - angular_min) / angular_step * (-maxY)
  207.         scale    = exp( log(rad_step)/rad_step ) ** (maxX)
  208.         scale*=prescale
  209.         rotation+=prerotate
  210.        
  211.         //result("\n\n\t Found shifts :x="+maxX+" y="+maxY)
  212.         //result("\n\t Found scale: "+scale+" ("+prescale+")")
  213.         //result("\n\t Found rotation : "+rotation/Pi()*180+" degree, or "+(rotation/Pi()*180-180)+" degree ("+preRotate/Pi()*180+")\n\n")
  214.     }
  215.    
  216.     void FindTranslationScaleAndRotation( object self, image ref, image sense, number &scale, number &rotation, number &ShiftX, number &ShiftY )
  217.     {
  218.         // First determine scale & rotation in a translational-invariate manner
  219.         result("\nInitial scale    : "+scale+"\n")
  220.         result("Initial rotation : "+rotation/Pi()*180+" [degree]\n")
  221.         result("Initial shifts XY: "+ShiftX+"/"+ShiftY+" [pixel]\n\n")
  222.        
  223.         OpenAndSetProgressWindow( "Find T&S&R","get rought S&R","")
  224.         self.FindScaleAndRotation( ref, sense, scale, rotation, shiftX, shiftY );
  225.        
  226.         OpenAndSetProgressWindow( "Find T&S&R","find T","S="+Format(scale,"%3.2f")+" R="+Format( rotation*180/Pi(), "%5.1f") )
  227.         image senseImage1 := self.ScaleRotateImage( sense, 1/scale, rotation )
  228.         senseImage1 := self.ZeroPadImage( ref, senseImage1, shiftX, shiftY )
  229.            
  230.         image senseImage2 := self.ScaleRotateImage( sense, 1/scale, rotation-Pi() )
  231.         senseImage2 := self.ZeroPadImage( ref, senseImage2, shiftX, shiftY)
  232.  
  233.         number maxX1, maxY1, maxV1
  234.         self.FindPhaseCorrelationMaximum( ref, senseImage1, maxX1, maxY1, maxV1, TRUE )
  235.         number maxX2, maxY2, maxV2
  236.         self.FindPhaseCorrelationMaximum( ref, senseImage2, maxX2, maxY2, maxV2, TRUE )
  237.         ShiftX += (maxV1>maxV2) ? maxX1 : maxX2
  238.         ShiftY += (maxV1>maxV2) ? maxY1 : maxY2
  239.         rotation = (maxV1>maxV2) ? rotation : rotation-Pi()
  240.         rotation += (rotation<0) ? 2*Pi() : 0
  241.         result("Found scale      : "+scale+"\n")
  242.         result("Found rotation   : "+rotation/Pi()*180+" [degree]\n")
  243.         result("Found shifts XY  : "+ShiftX+"/"+ShiftY+" [pixel]\n")
  244.         result("Phase correlation: "+max(maxV1,maxV2)+"\n")
  245.         OpenAndSetProgressWindow( "","","")
  246.     }
  247.    
  248.     ///////////////////////
  249.     // Results & Display //
  250.     ///////////////////////
  251.     image RGBOverlay( object self, image ref, image sense, number scale, number rotation, number shftX, number shftY, number &getCorMax )
  252.     {
  253.         image sen := self.ScaleRotateImage( sense, 1/scale, rotation )
  254.         image cut := self.ZeroPadImage( ref, sen, shftX, shftY )
  255.         if ( TRUE == getCorMax )
  256.         {
  257.             image ref_cut = ( cut == 0 ) ? 0 : ref
  258.             image crossCorr := CrossCorrelate( cut, ref_cut )
  259.             getCorMax = max(crossCorr)
  260.             Result( "\n Correlation: " + getCorMax  +"\n")
  261.         }
  262.            
  263.         return RGB(ref/max(ref)*255,cut/max(cut)*255,0)
  264.     }
  265.    
  266.     number GetCorrelationCoefficient( object self, image ref, image sense,  number scale, number rotation, number shftX, number shftY )
  267.     {
  268.         image sen := self.ScaleRotateImage( sense, 1/scale, rotation )
  269.         image sen_cut := self.ZeroPadImage( ref, sen, shftX, shftY )
  270.         image ref_cut = ( sen_cut == 0 ) ? 0 : ref
  271.         image crossCorr := CrossCorrelate( sen_cut, ref_cut )
  272.         return max(crossCorr)
  273.     }
  274. }
  275.  
  276.  
  277. //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  278. // MAIN PROGRAM (Testing RegisterClass)                                                                                                                         //
  279. //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  280. {
  281.     kDebug = 0
  282.     image ref := B
  283.     image sen := D
  284.     object ob = Alloc(RegisterSubImage)
  285.     number mx,my,mv
  286.     number scal,rot,shftx,shfty
  287.     // If a rough-aligment has already been made,
  288.     // it can be fed in here to improve results of the algorithm.
  289.     scal    = 4//1.7 //1.7
  290.     rot     = 0//-45 *PI()/180
  291.     shftx   = 0//201
  292.     shfty   = 0//270
  293.     //ob.ScaleRotateImage( sen, 1/scal, rot ).showimage()
  294.     //exit(0)
  295.     number CorMax = TRUE
  296.     ClearResults()
  297.     ob.FindTranslationScaleAndRotation(ref,sen,scal,rot,shftx,shfty )
  298.     ob.RGBOverlay(ref,sen,scal,rot,shftx,shfty,CorMax).ShowImage()
  299. }
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×