[AstroPy] Astropy for imaging

[R Schumacher]

At 02:53 PM 5/14/2015, you wrote:
>  I hate to ask too much but I haven't been able to find good material for c
>  oming up with image stacking algorithms of my own, would anybody 
> be able to point me to some material that covers that in detail, 
> preferably with Python examples?

Take a look at  Xie's Matlab work, which I translated to Python
http://rjs.org/astro/1004x/Python/register/
and another paper for log-polar registration (ie scale/shift)
https://etd.ohiolink.edu/rws_etd/document/get/osu1236610454/inline

Jan Kybic
http://cmp.felk.cvut.cz/~kybic/
has also published advanced registration methods, some very recently. 
His older "image registration by warping" code is archived at
https://web.archive.org/web/20041209143140/http://cmp.felk.cvut.cz/~kybic/thesis/jkpyreg-0.2.tgz

https://github.com/KitwareMedical/pyLAR

http://www0.cs.ucl.ac.uk/opensource_mia_ws_2012/links.html#registration



Ray Schumacher
Programmer/Consultant
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20150515/010527db/attachment.html>
-------------- next part --------------
"""
;  Authors:
;  Hongjie Xie,  Nigel Hicks,  G. Randy Keller, Haitao Huang, Vladik Kreinovich
;
;  Title of the paper:
;  An IDL/ENVI implementation of the FFT Based Algorithm for Automatic Image Registration


;  IDL source code for shift only
;  ----------------------------------------------------------------------------
"""


import numarray
import numarray.fft as fft
#import numarray.nd_image as nd_image
#import math
import string
import sys

def SizeImage(FileName):
    hdr_lun1 = open(FileName+'.hdr', 'r')
    info = hdr_lun1.readlines()
    hdr_lun1.close()
    samples = int(string.split(info[3])[2])
    lines = int(string.split(info[4])[2])
    return (samples, lines)

def LoadImage(FileName, samples, lines):
    ## cols=samples, rows=lines
    data_lun1 = open(FileName+'.img', 'rb')
    IArray = numarray.fromstring(data_lun1.read(), numarray.UInt8, shape=(lines, samples)).astype(numarray.Float)
    data_lun1.close()

    return numarray.transpose(IArray)


##=====================================================================
##
## MAIN
##
##=====================================================================

print 'Loading images...'
FileName = 't350380ori'
#FileName = '6cx7r'
cols, rows = SizeImage(FileName)
FFTnorm = complex(cols*rows, 0)
print 'cols=', cols ,'rows=',rows
Im1 = LoadImage(FileName, cols, rows)

FileName = 't350380shf'
#FileName = '6cx7r'
cols, rows = SizeImage(FileName)
Im2 = LoadImage(FileName, cols, rows)
print Im1[0,0], Im1[1,0], Im1[0,1], Im1[cols/2,rows/2]
print ''

print 'Doing FFT on two images...'
Ffft_im1=fft.fft2d(Im1)/FFTnorm
Ffft_im2=fft.fft2d(Im2)/FFTnorm
print Ffft_im1[0,0], Ffft_im1[cols/2,rows/2]

print 'doing the ratio'
R1= (Ffft_im1 * Ffft_im2.conjugate())
print R1[0,0], R1[1,0], R1[cols/2,rows/2]
R2= (abs(Ffft_im1) * abs(Ffft_im2))
print R2[0,0], R2[cols/2,rows/2]
R = R1 / R2
print R[0,0], R[cols/2,rows/2]

print 'doing inverse_fft2d'
IR=fft.inverse_fft2d(R) * FFTnorm
print IR[0,0], IR[cols/2,rows/2]

flat_IR = numarray.ravel(numarray.transpose(IR))
print flat_IR.type()
I= numarray.argmax(flat_IR)
maxn=flat_IR[I]


IX= I % cols
IY= I / cols

X=IX
Y=IY
print 'I=', I, 'X=', X, 'Y=', Y


if((X > 0.5*cols) and ( Y > 0.5*rows)):
    X=X-cols
    Y=Y-rows
    IX=X
    IY=Y
elif(X > 0.5*cols):
    X=X-cols
    IX=X
    IY=IY
elif(Y > 0.5*rows):
    Y=Y-rows
    IX=IX
    IY=Y
else:
    IX=IX
    IY=IY


print 'this is max:', maxn
print 'position: x=', IX, ', y=', IY