source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_TubeAdjustments.ipf @ 1037

Last change on this file since 1037 was 1024, checked in by srkline, 6 years ago

minor changes to prefix functions with "V_" to avoid conflicts with non-VSANS functions.

File size: 13.7 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3//
4// functions for testing and then actually applying the nonlinear corrections to the
5// tube detectors. These routines are for a test bank of 8 tubes (vertical) that were
6// run at a subdivision of 1024. VSANS will be different in practice
7//
8// but the fundamental process is the same, and can be translated into proper functions as needed
9//
10
11
12
13// the main routines are:
14
15//(1)
16//to get from individual tubes to an array
17//      Tubes_to_Array()                       
18
19//(2)
20// then to locate all of the peak positions
21//      MakeTableForPeaks(numTube,numPeak)             
22//      Identify_AllPeaks()
23//              AutoFindPeaksCustom()           // if Identify_AllPeaks  doesn't work -try this, setting the "noise" to 1 and smoothing to 2
24
25//(3)
26// fit to find all of the quadratic coefficients
27//      MakeTableForFitCoefs(numTube,numCoef)
28//      PlotFit_AllPeaks()
29
30
31//(4)
32// then pick a display method
33//
34//      Make_mm_tubes()
35//      Append_Adjusted_mm()
36//
37//      MakeMatrix_PixCenter()
38//      FillMatrix_Pix_Center(ind)
39//
40//
41// -or- (note that the pack_image wave that is generated here is for display ONLY)
42// --since it is interpolated
43//
44// Interpolate_mm_tubes()
45
46
47// The function most used externally is:
48// V_TubePix_to_mm(c0,c1,c2,pix)
49//
50// which will return the real space distance (in mm?) for a given pixel
51// and the tube's coefficients. The distance is relative to the zero position on the
52// detector (which was defined when the coefficients were determined)
53
54
55
56
57
58
59
60// (0) -- what I start with:
61// -- a table of the mm spacing of the slots (20 of them)
62// -- masked data from each of the (8) tubes
63// -- the table of slots may need to be corrected for parallax, depending on the geometry of the test
64// ** In the table of slots, pick a slot near the center, and SET that to ZERO. Then all of the other
65//   distances are relative to that zero point. This is a necessary reference point.
66//
67
68
69// (1) -- get the individual tubes into an array
70//
71//
72Proc V_Tubes_to_Array()
73        Make/O/D/N=(8,1127) pack
74        edit pack
75        display;appendimage pack
76        pack[0][] = tube1[q]
77        pack[1][] = tube2[q]
78        pack[2][] = tube3[q]
79        pack[3][] = tube4[q]
80        pack[4][] = tube5[q]
81        pack[5][] = tube6[q]
82        pack[6][] = tube7[q]
83        pack[7][] = tube8[q]
84        ModifyImage pack ctab= {*,*,ColdWarm,0}
85End
86
87// (2) -- for each of the tubes, find the x-position (in pixels) of each of the (20) peaks
88// -- load the Analysis Package "MultiPeakFit 2"
89//
90// automatically find the peaks (after including MultiPeakFit 2)
91//              AutomaticallyFindPeaks()
92//
93//-- or if having difficulty
94//              AutoFindPeaksCustom()           // try this, setting the "noise" to 1 and smoothing to 2
95//
96// -- if really having difficulty, you can do the "full" MultiPeak Fit
97//
98// -- If (hopefully) using the easy way, the results are in:
99// root:WA_PeakCentersY,root:WA_PeakCentersX
100//
101// -- so to see the results:
102//¥Edit/K=0  root:WA_PeakCentersY,root:WA_PeakCentersX
103//
104// -- then sort the results - they seem to be in no real order...
105//¥Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
106//
107Proc V_MakeTableForPeaks(numTube,numPeak)
108        Variable numTube,numPeak
109       
110        Make/O/D/N=(numPeak,numTube) PeakTableX,peakTableY              //*2 to store x-location and peak height (y)
111        Edit peakTableX
112End
113
114Proc V_Identify_AllPeaks()
115
116        Variable ii,numTubes=8
117        String str="tube"
118       
119        ii=1
120        do
121                V_Identify_Peaks(str+num2str(ii),ii-1)
122                ii+=1
123        while(ii<=numTubes)
124
125End
126
127Proc V_Identify_Peaks(tubeStr,ind)
128        String tubeStr
129        Variable ind
130       
131        // must use a wave of pixels rather than "calculated" -- if calculated is used it only
132        // returns integer values for the peak locations
133       
134//      AutomaticallyFindPeaks() //-- this is a function that doesn't take any parameters - so
135// I need to pull this from the WM function to call the worker directly
136        Variable pBegin=0, pEnd= numpnts($(tubeStr))-1
137        Variable/C estimates= EstPeakNoiseAndSmfact($(tubeStr),pBegin, pEnd)
138        Variable noiselevel=real(estimates)
139        Variable smoothingFactor=imag(estimates)
140        Variable maxPeaks = 20
141        Variable minPeakPercent = 10
142       
143        AutoFindPeaksWorker($(tubeStr), $("tube_pixel"), pBegin, pEnd, maxPeaks, minPeakPercent, noiseLevel, smoothingFactor)
144// end WM function call
145
146        Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
147       
148        peakTableX[][ind] = WA_PeakCentersX[p]
149        peakTableY[][ind] = WA_PeakCentersY[p]
150       
151End
152
153// -- save a copy of the root:WA_PeakCentersY,root:WA_PeakCentersX values
154//    for later in case the fitting failed, then you can go back and re-do
155//
156// -- then plot:
157//
158//      Display peak_spacing_mm_ctr vs WA_PeakCentersX
159//
160//      Then do a "QuickFit" of the peak position to the data using a polynomial of order 3 (= quadratic)
161//
162// result is in W_coef, W_sigma
163//
164// -- an example of the "quickFit" command is below, so it can be programmed rather than the menu every time
165//¥Display peak_spacing_mm_ctr vs WA_PeakCentersX3
166//¥CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=WA_PeakCentersX3/D
167//  fit_peak_spacing_mm_ctr= poly(W_coef,x)
168//  W_coef={-571.42,1.1135,-4.2444e-05}
169//  V_chisq= 8.5841;V_npnts= 20;V_numNaNs= 0;V_numINFs= 0;
170//  V_startRow= 0;V_endRow= 19;
171//  W_sigma={0.595,0.00246,2.15e-06}
172//  Coefficient values ± one standard deviation
173//      K0      =-571.42 ± 0.595
174//      K1      =1.1135 ± 0.00246
175//      K2      =-4.2444e-05 ± 2.15e-06
176//
177//
178//
179// for (8) tubes, keep all of the fit coefficients
180//
181//¥make/O/D/N=(3,8) fit_coef
182//¥edit fit_coef
183//¥make/O/D/N=(3,8) fit_sigma
184//¥edit fit_sigma
185//
186// -- copy and paste in the W_coef and W_sigma values (or by a command)
187//
188
189
190Proc V_MakeTableForFitCoefs(numTube,numCoef)
191        Variable numTube,numCoef
192       
193        Make/O/D/N=(numCoef,numTube) TubeCoefTable,TubeSigmaTable               //
194        Edit TubeCoefTable
195End
196
197Proc V_PlotFit_AllPeaks()
198
199        Variable ii,numTubes=8
200       
201        ii=1
202        do
203                V_PlotFit_Peaks(ii-1)
204                ii+=1
205        while(ii<=numTubes)
206
207End
208
209Proc V_PlotFit_Peaks(ind)
210        Variable ind
211       
212        //hopefully 20 points - need better control of this
213        Duplicate/O WA_PeakCentersX, tmpX
214       
215        tmpX = peakTableX[p][ind]
216        Display peak_spacing_mm_ctr vs tmpX
217       
218        CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=tmpX/D
219       
220        TubeCoefTable[][ind] = W_coef[p]
221        TubeSigmaTable[][ind] = W_sigma[p]
222       
223End
224
225
226
227
228
229//¥Duplicate tube1 tube1_mm
230//¥tube1_mm = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],p)
231
232
233
234
235////////
236// then there are various display options:
237
238// adjust the center (pixel) of the tube:
239// - measCtr is the pixel location of the DEFINED "zero" peak
240// nominal Ctr is the pixel number of this DEFINED zer0 position, nominally nPix/2-1
241//
242// ( be sure to pick better names, and use a loop...)
243//      adj_tube = raw_tube[p+(measCtr-nominalCtr)]
244//
245// then fill and display a new matrix. The center will be reasonably well aligned, and will
246// get worse towards the ends of the tubes
247// (this may be the "preferred" way of displaying the raw data if the centers are far off)
248// -- this also may be what I need to start with to fit the data to locate the beam center.
249
250
251// I can also display the fully corrected tubes, where the y-axis is now in real-space mm
252// rather than arbitrary pixels. The x-axis is still tube nubmer.
253// -- do this with the procedure"
254//   Append_Adjusted_mm()               // name may change...
255//
256// -- the point of the appending is that it allows each tube to be plotted on an image plot
257// with its own y-axis. Every one will be different and will be non-linear. These conditions
258// BOTH prevent using any of the "normal" image plotting or manipulation routines.
259// - the gist is below:
260//
261//      duplicate tube1_mm adjusted_mm_edge
262//      InsertPoints 0,1, adjusted_mm_edge
263//      // be sure to use the correct set of coefficients
264//      adjusted_mm_edge[0] = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],-1)
265//     
266//      Display;AppendImage adjusted_pack vs {*,adjusted_mm_edge}
267
268
269Proc V_Make_mm_tubes()
270
271        Variable ii,numTubes=8
272       
273        ii=1
274        do
275                Duplicate $("tube"+num2str(ii)) $("tube"+num2str(ii)+"_mm")
276                $("tube"+num2str(ii)+"_mm") = V_TubePix_to_mm(TubeCoefTable[0][ii-1],TubeCoefTable[1][ii-1],TubeCoefTable[2][ii-1],p)
277                ii+=1
278        while(ii<=numTubes)
279       
280End
281
282
283Proc V_Append_Adjusted_mm()
284
285// a blank base image
286        Duplicate/O pack junk
287        junk=0
288        SetScale/I y -600,600,"", junk          // -600,600 is tooo large and not general
289        Display;Appendimage junk
290
291        Variable ii
292       
293        ii=1
294        do
295                make/O/D/N=(1,1127) $("tube"+num2str(ii)+"_mm_mat")=0   
296       
297                $("tube"+num2str(ii)+"_mm_mat")[0][] = $("tube"+num2str(ii))
298                SetScale/I x (ii-1),ii,"", $("tube"+num2str(ii)+"_mm_mat")              //builds up the x-axis
299               
300                duplicate/O $("tube"+num2str(ii)+"_mm") $("edge"+num2str(ii)+"_mm")
301                InsertPoints 0,1, $("edge"+num2str(ii)+"_mm")           //needs to be one point longer
302        // be sure to use the correct set of coefficients
303                $("edge"+num2str(ii)+"_mm")[0] = V_TubePix_to_mm(TubeCoefTable[0][0],TubeCoefTable[1][0],TubeCoefTable[2][0],-1)
304       
305                AppendImage $("tube"+num2str(ii)+"_mm_mat") vs {*,$("edge"+num2str(ii)+"_mm")}
306                ModifyImage $("tube"+num2str(ii)+"_mm_mat") ctab= {*,*,ColdWarm,0}
307       
308                ii+=1
309        while(ii < 9)
310       
311end
312
313
314
315////////////////////////////////
316
317Proc V_MakeMatrix_PixCenter()
318        Duplicate/O pack pack_centered
319       
320        Variable ii,numTubes=8
321       
322        ii=1
323        do
324                V_FillMatrix_Pix_Center(ii)
325                ii+=1
326        while(ii<=numTubes)
327       
328        Display;AppendImage pack_centered
329        ModifyImage pack_centered ctab= {*,*,ColdWarm,0}
330
331end
332
333//
334// this fills a matrix with the tubes center aligned ONLY, with the y-axis in pixels
335//
336// adj_tube = raw_tube[p+(measCtr-nominalCtr)]
337// finds the center automatically from the tubeN_mm wave, where it crosses zero
338//
339// Tube #1 is the "base" ans others are shifted to match that one's "zero"
340//
341// FindRoots/P=W_coef           can also be used to find the roots (but which one?)
342//
343Function V_FillMatrix_Pix_Center(ind)
344        Variable ind
345       
346        Wave adj=root:pack_centered
347        Wave tube1_mm = $("root:tube1_mm")
348        Wave tube = $("root:tube"+num2str(ind))
349        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
350
351        Variable base,shift,ii,num,pt
352       
353        num=numpnts(tube)
354       
355        FindLevel tube1_mm 0
356        base = round(V_LevelX)
357       
358       
359        FindLevel tube_mm 0
360        shift = round(V_LevelX)
361       
362        for(ii=0;ii<num;ii+=1)
363                pt = ii + (shift-base)
364                if(pt >= 0 && pt < num)
365                        adj[ind-1][ii] = tube[pt]
366                endif
367        endfor
368       
369        return(0)
370End
371
372
373// this fills a matrix with the tubes center aligned ONLY, with the y-axis in mm
374// -- there seems to be little reason to do this --
375// -- either keep pixels and align centers
376// -- OR -- use mm and append each tube with its own y-axis
377//
378Function V_FillAdjusted(ind)
379        Variable ind
380       
381        Wave adj=root:adjusted_pack
382        Wave tube1_mm
383        Wave tube = $("root:tube"+num2str(ind))
384        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
385
386        Variable base,shift,ii,num,pt
387       
388        num=numpnts(tube1_mm)
389       
390        FindLevel tube1_mm 0
391        base = round(V_LevelX)
392       
393       
394        FindLevel tube_mm 0
395        shift = round(V_LevelX)
396       
397        for(ii=0;ii<num;ii+=1)
398                pt = ii + (shift-base)
399                if(pt >= 0 && pt < num)
400                        adj[ind-1][ii] = tube[pt]
401                endif
402        endfor
403       
404        return(0)
405End
406
407
408// c0,c1,c2,pix
409// c0-c2 are the fit coefficients
410// pix is the test pixel
411//
412// returns the distance in mm (relative to ctr pixel)
413// ctr is the center pixel, as defined when fitting to quadratic
414//
415Function V_TubePix_to_mm(c0,c1,c2,pix)
416        Variable c0,c1,c2,pix
417       
418        Variable dist
419        dist = c0 + c1*pix + c2*pix*pix
420       
421        return(dist)
422End
423
424////
425
426
427
428// set the (linear) range of the y-axis of the matrix to be the
429// range of the 1st tube. This is completely arbitrary
430//
431Proc V_Interpolate_mm_tubes()
432
433        Duplicate/O pack pack_image
434
435        Variable ii,numTubes=8
436        Variable p1,p2
437        p1 = tube1_mm[0]
438        p2 = tube1_mm[numpnts(tube1_mm)-1]
439       
440        SetScale/I y p1,p2,"", pack_image
441       
442        // then make a temporary 1D wave to help with the interpolation
443        Duplicate/O tube1_mm lin_mm,lin_val
444        SetScale/I x p1,p2,"", lin_mm
445        lin_mm = x                      //fill with the linear mm spacing
446        lin_val=0
447       
448        ii=1
449        do
450                lin_val = interp(lin_mm, $("tube"+num2str(ii)+"_mm"), $("tube"+num2str(ii)))
451                pack_image[ii-1][] = lin_val[q]
452               
453                ii+=1
454        while(ii<=numTubes)
455       
456        display;appendimage pack_image
457        ModifyImage pack_image ctab= {*,*,ColdWarm,0}
458       
459End
460
461
462
463
464
465
466
467
468
469
470
471
472
473// this doesn't work - the interpolation step doesn't do what I want
474// and the plot of the triplet with f(z) for color doesn't fill space like I want
475Proc V_AnotherExample()
476
477        Concatenate/O/NP {tube1_mm,tube2_mm,tube3_mm,tube4_mm,tube5_mm,tube6_mm,tube7_mm,tube8_mm},cat_mm
478        Concatenate/O/NP {tube1,tube2,tube3,tube4,tube5,tube6,tube7,tube8},cat_tubes
479        Duplicate/O cat_mm,cat_num
480        Variable num=1127
481        cat_num[0,num-1]=1
482        cat_num[num,2*num-1]=2
483        cat_num[2*num,3*num-1]=3
484        cat_num[3*num,4*num-1]=4
485        cat_num[4*num,5*num-1]=5
486        cat_num[5*num,6*num-1]=6
487        cat_num[6*num,7*num-1]=7
488        cat_num[7*num,8*num-1]=8
489
490        Display cat_mm vs cat_num
491        ModifyGraph mode=3,marker=9
492        ModifyGraph zColor(cat_mm)={cat_tubes,*,*,ColdWarm,0}
493
494        Concatenate/O {cat_num,cat_mm,cat_tubes}, tripletWave
495        ImageInterpolate Kriging tripletWave
496        AppendImage M_InterpolatedImage
497
498//      Make/O/N=20 xWave=enoise(4),yWave=enoise(5),zWave=enoise(6)  // Random points
499//      Display yWave vs xWave
500//      ModifyGraph mode=3,marker=19
501//      ModifyGraph zColor(yWave)={zWave,*,*,Rainbow,0}
502//
503//      Concatenate/O {xWave,yWave,zWave}, tripletWave
504//      ImageInterpolate/S={-5,0.1,5,-5,0.1,5} voronoi tripletWave
505//      AppendImage M_InterpolatedImage
506
507end
508
509// this desn't work either...
510// (same y-axis for the entire image, which is not the case for the tubes)
511//
512// from the WM help file:
513// Plotting a 2D Z Wave With 1D X and Y Center Data
514//
515Function V_MakeEdgesWave(centers, edgesWave)
516        Wave centers    // Input
517        Wave edgesWave  // Receives output
518       
519        Variable N=numpnts(centers)
520        Redimension/N=(N+1) edgesWave
521
522        edgesWave[0]=centers[0]-0.5*(centers[1]-centers[0])
523        edgesWave[N]=centers[N-1]+0.5*(centers[N-1]-centers[N-2])
524        edgesWave[1,N-1]=centers[p]-0.5*(centers[p]-centers[p-1])
525End
526
527//This function demonstrates the use of MakeEdgesWave:
528Function V_DemoPlotXYZAsImage()
529        Make/O mat={{0,1,2},{2,3,4},{3,4,5}}    // Matrix containing Z values
530        Make/O centersX = {1, 2.5, 5}           // X centers wave
531        Make/O centersY = {300, 400, 600}               // Y centers wave
532        Make/O edgesX; V_MakeEdgesWave(centersX, edgesX)        // Create X edges wave
533        Make/O edgesY; V_MakeEdgesWave(centersY, edgesY)        // Create Y edges wave
534        Display; AppendImage mat vs {edgesX,edgesY}
535End
536
Note: See TracBrowser for help on using the repository browser.