Hi Sepp, I did do a calibration script on a micro CT platform. The idea was that I had a bb taken from 4 different view points at 4 different heights with a known distance. Their projected position was segmented beforehand on the projections and I tried to match the projections by optimizing the calibration (which we see in before.pdf and after.pdf). I don't think it's a good script, e.g., because I tried to find 8 parameters when 7 are sufficient for my system (see, e.g., this paper <http://iopscience.iop.org/article/10.1088/0031-9155/45/11/327/meta;jsessionid=8DD0F1160098FBC65BB161A0257FBA75.c3.iopscience.cld.iop.org#>). But this gives you a starting point to develop your own script which will hopefully be better. Please share if you do something from it! Thanks, Simon
On Sat, Nov 19, 2016 at 8:51 PM, Sepp de Raedt <[email protected]> wrote: > Hi RTK users, > > I'm interested in reconstructing images acquired by a fluoroscopy system > rotating around a knee phantom. Unfortunately, I've been unsuccessful so > far. I recognize some of the phantom in the reconstructed image, but it > contains double contours and the shape is deformed. > > I think the issue might be incorrect geometry specification. I have images > also including a phantom containing markers for which we know the > coordinates, which should allow us to calibrate the system. On the wiki on > image quality, I saw that some users had developed scripts to do the > calibration? Would someone be willing to share it? > > Kind regards, > Sepp > > ______________________________________________________________________ > This email has been scanned by the Symantec Email Security.cloud service. > For more information please visit http://www.symanteccloud.com > ______________________________________________________________________ > _______________________________________________ > Rtk-users mailing list > [email protected] > http://public.kitware.com/mailman/listinfo/rtk-users >
after.pdf
Description: Adobe PDF document
before.pdf
Description: Adobe PDF document
#!/usr/bin/env python
from __future__ import print_function
import SimpleITK as sitk
import SimpleRTK as srtk
import sys
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
from scipy.optimize import minimize
import argparse
import csv
import scipy.optimize
segpos = np.loadtxt('CentroidPositions_H5_to_H25_90degstep.csv',delimiter=',')
# Convert pixel positions to mm positions
segpos[:,2] = (segpos[:,2]-1023.5)*0.01709
segpos[:,3] = (segpos[:,3]-1023.5)*0.01709
# Parameters: angle1,angle2,bbx,bbz,sid,sdd,px,py,ia,oa,sx,sy
optp = np.array([0,0,0,0,0,0,0,0,0,0,0,0])
def GetProjPos(optp):
geo = srtk.ThreeDCircularProjectionGeometry()
resultx = np.zeros(segpos.shape[0])
resulty = np.zeros(segpos.shape[0])
for i in range(4):
geo.AddProjection(210+optp[4],270+optp[5],-90.*i,optp[6],optp[7],optp[8],optp[9],optp[10],optp[11])
for i in range(segpos.shape[0]):
igeo = np.round(segpos[i,0]/90).astype('int')
pos3D = np.zeros(4)
pos3D[0] = optp[0] + (-15+segpos[i,1])*np.sin(optp[2]*np.pi/180.)*np.sin(optp[3]*np.pi/180.)
pos3D[1] = (-15+segpos[i,1])*np.cos(optp[2]*np.pi/180.)
pos3D[2] = optp[1] + (-15+segpos[i,1])*np.sin(optp[2]*np.pi/180.)*np.cos(optp[3]*np.pi/180.)
pos3D[3] = 1.
projMat = np.reshape(np.array(geo.GetMatrix(igeo)), [3,4])
projPos = np.matmul(projMat,pos3D)
projPos /= projPos[2]
resultx[i] = projPos[0]
resulty[i] = projPos[1]
#return np.concatenate((resultx,resulty))
return segpos[:,2],resultx,segpos[:,3],resulty
def Plot(optp, name):
plt.clf()
sx,rx,sy,ry = GetProjPos(optp)
plt.plot(sx,sy, 'ko', markersize=5)
plt.plot(rx,ry, 'ro', markersize=5)
plt.savefig(name, bbox_inches='tight')
def CostFunc(optp):
sx,rx,sy,ry = GetProjPos(optp)
return np.sum((sx-rx)**2+(sy-ry)**2)
#print(scipy.optimize.leastsq(CostFunc, optp, full_output = 0, args=()))
res = minimize(CostFunc, optp, args=(), method='Powell', options={'disp':True,'maxiter':1e6,'maxfev': 1e6, 'xtol': 1e-15, 'ftol': 1e-15})
Plot(optp, 'before.pdf')
Plot(res.x, 'after.pdf')
gfull = srtk.ThreeDCircularProjectionGeometry()
for i in range(360):
gfull.AddProjection(210+res.x[4],270+res.x[5],-1.*i,res.x[6],res.x[7],res.x[8],res.x[9],res.x[10],res.x[11])
geometryWriter = srtk.ThreeDCircularProjectionGeometryXMLFileWriter()
geometryWriter.SetFileName('gfull')
geometryWriter.Execute(gfull)
_______________________________________________ Rtk-users mailing list [email protected] http://public.kitware.com/mailman/listinfo/rtk-users
