43 from shesha.util import dm_util, influ_util, kl_util
49 from scipy
import interpolate
52 from typing
import List
57 def dm_init(context: carmaWrap_context, p_dms: List[conf.Param_dm],
58 p_tel: conf.Param_tel, p_geom: conf.Param_geom,
59 p_wfss: List[conf.Param_wfs] =
None, keepAllActu: bool =
False) -> Dms:
60 """Create and initialize a Dms object on the gpu
63 context: (carmaWrap_context): context
64 p_dms: (list of Param_dms) : dms settings
65 p_tel: (Param_tel) : telescope settings
66 p_geom: (Param_geom) : geom settings
67 p_wfss: (list of Param_wfs) : wfs settings
69 Dms: (Dms): Dms object
72 if (p_wfss
is not None):
75 for i
in range(len(p_wfss)):
76 xpos_wfs.append(p_wfss[i].xpos)
77 ypos_wfs.append(p_wfss[i].ypos)
83 types_dm = [p_dm.type
for p_dm
in p_dms]
84 if scons.DmType.TT
in types_dm:
85 first_TT = types_dm.index(scons.DmType.TT)
86 if np.any(np.array(types_dm[first_TT:]) != scons.DmType.TT):
87 raise RuntimeError(
"TT must be defined at the end of the dms parameters")
89 for i
in range(len(p_dms)):
90 max_extent = _dm_init(context, dms, p_dms[i], xpos_wfs, ypos_wfs, p_geom,
91 p_tel.diam, p_tel.cobs, p_tel.pupangle, max_extent,
92 keepAllActu=keepAllActu)
97 def _dm_init(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm, xpos_wfs: list,
98 ypos_wfs: list, p_geom: conf.Param_geom, diam: float, cobs: float,
99 pupAngle: float, max_extent: int, keepAllActu: bool =
False):
100 """ inits a Dms object on the gpu
103 context: (carmaWrap_context): context
104 dms: (Dms) : dm object
106 p_dm: (Param_dms) : dm settings
108 xpos_wfs: (list) : list of wfs xpos
110 ypos_wfs: (list) : list of wfs ypos
112 p_geom: (Param_geom) : geom settings
114 diam: (float) : diameter of telescope
116 cobs: (float) : cobs of telescope
118 max_extent: (int) : maximum dimension of all dms
121 max_extent: (int) : new maximum dimension of all dms
125 if (p_dm.pupoffset
is not None):
126 p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
128 patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
131 if (p_dm.type == scons.DmType.PZT):
132 if p_dm.file_influ_hdf5 ==
None:
133 p_dm._pitch = patchDiam / float(p_dm.nact - 1)
135 extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
136 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent,
140 make_pzt_dm(p_dm, p_geom, cobs, pupAngle, keepAllActu=keepAllActu)
145 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
147 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
148 p_dm._dim_screen = dim
149 ninflupos = p_dm._influpos.size
150 n_npts = p_dm._ninflu.size
151 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
152 ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
154 dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
155 p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
157 elif (p_dm.type == scons.DmType.TT):
159 if (p_dm.alt == 0)
and (max_extent != 0):
160 extent = int(max_extent * 1.05)
161 if (extent % 2 != 0):
164 extent = p_geom.pupdiam + 16
165 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
167 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
169 dim = p_dm._n2 - p_dm._n1 + 1
171 p_dm._dim_screen = dim
172 dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
173 context.active_device)
174 dms.d_dms[-1].tt_loadarrays(p_dm._influ)
176 elif (p_dm.type == scons.DmType.KL):
178 extent = p_geom.pupdiam + 16
179 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
181 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
183 dim = p_dm._n2 - p_dm._n1 + 1
188 p_dm._dim_screen = dim
190 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
191 p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
193 dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
198 raise TypeError(
"This type of DM doesn't exist ")
206 def _dm_init_factorized(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm,
207 xpos_wfs: list, ypos_wfs: list, p_geom: conf.Param_geom,
208 diam: float, cobs: float, pupAngle: float, max_extent: int,
209 keepAllActu: bool =
False):
210 """ inits a Dms object on the gpu
214 context: (carmaWrap_context): context
215 dms: (Dms) : dm object
217 p_dm: (Param_dms) : dm settings
219 xpos_wfs: (list) : list of wfs xpos
221 ypos_wfs: (list) : list of wfs ypos
223 p_geom: (Param_geom) : geom settings
225 diam: (float) : diameter of telescope
227 cobs: (float) : cobs of telescope
229 max_extent: (int) : maximum dimension of all dms
232 max_extent: (int) : new maximum dimension of all dms
236 if (p_dm.pupoffset
is not None):
237 p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
239 patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
242 if (p_dm.type == scons.DmType.PZT)
and p_dm.file_influ_hdf5
is not None:
245 if (p_dm.type == scons.DmType.PZT):
246 p_dm._pitch = patchDiam / float(p_dm.nact - 1)
248 extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
251 make_pzt_dm(p_dm, p_geom, cobs, pupAngle, keepAllActu=keepAllActu)
253 elif (p_dm.type == scons.DmType.TT):
254 if (p_dm.alt == 0)
and (max_extent != 0):
255 extent = int(max_extent * 1.05)
256 if (extent % 2 != 0):
259 extent = p_geom.pupdiam + 16
261 elif (p_dm.type == scons.DmType.KL):
262 extent = p_geom.pupdiam + 16
264 raise TypeError(
"This type of DM doesn't exist ")
270 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
273 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
275 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
277 if (p_dm.type == scons.DmType.PZT):
278 ninflupos = p_dm._influpos.size
279 n_npts = p_dm._ninflu.size
280 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
281 ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
283 dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
284 p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
285 elif (p_dm.type == scons.DmType.TT):
287 dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
288 context.active_device)
289 dms.d_dms[-1].tt_loadarrays(p_dm._influ)
290 elif (p_dm.type == scons.DmType.KL):
294 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
295 p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
296 dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
302 def dm_init_standalone(context: carmaWrap_context, p_dms: list, p_geom: conf.Param_geom,
303 diam=1., cobs=0., pupAngle=0., wfs_xpos=[0], wfs_ypos=[0]):
304 """Create and initialize a Dms object on the gpu
307 p_dms: (list of Param_dms) : dms settings
309 p_geom: (Param_geom) : geom settings
311 diam: (float) : diameter of telescope (default 1.)
313 cobs: (float) : cobs of telescope (default 0.)
315 pupAngle: (float) : pupil rotation angle (degrees, default 0.)
317 wfs_xpos: (array) : guide star x position on sky (in arcsec).
319 wfs_ypos: (array) : guide star y position on sky (in arcsec).
323 if (len(p_dms) != 0):
325 for i
in range(len(p_dms)):
326 _dm_init(context, dms, p_dms[i], wfs_xpos, wfs_ypos, p_geom, diam, cobs,
327 pupAngle, max_extent)
331 def make_pzt_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, cobs: float,
332 pupAngle: float, keepAllActu: bool =
False):
333 """Compute the actuators positions and the influence functions for a pzt DM.
334 NOTE: if the DM is in altitude, central obstruction is forced to 0
337 p_dm: (Param_dm) : dm parameters
339 p_geom: (Param_geom) : geometry parameters
341 cobs: (float) : telescope central obstruction
344 influ: (np.ndarray(dims=3, dtype=np.float64)) : cube of the IF for each actuator
349 coupling = p_dm.coupling
356 if (p_dm.influ_type == scons.InfluType.PETAL):
360 if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
361 smallsize = influ_util.makeRadialSchwartz(pitch, coupling)
362 elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
363 smallsize = influ_util.makeSquareSchwartz(pitch, coupling)
364 elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
365 smallsize = influ_util.makeBlacknutt(pitch, coupling)
366 elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
367 smallsize = influ_util.makeGaussian(pitch, coupling)
368 elif (p_dm.influ_type == scons.InfluType.BESSEL):
369 smallsize = influ_util.makeBessel(pitch, coupling, p_dm.type_pattern)
370 elif (p_dm.influ_type == scons.InfluType.DEFAULT):
371 smallsize = influ_util.makeRigaut(pitch, coupling)
373 print(
"ERROR influtype not recognized ")
374 p_dm._influsize = smallsize
379 if p_dm.type_pattern
is None:
380 p_dm.type_pattern = scons.PatternType.SQUARE
382 if p_dm.type_pattern == scons.PatternType.HEXA:
383 print(
"Pattern type : hexa")
384 cub = dm_util.createHexaPattern(pitch, p_geom.pupdiam * 1.1)
386 elif p_dm.type_pattern == scons.PatternType.HEXAM4:
387 print(
"Pattern type : hexaM4")
389 cub = dm_util.createDoubleHexaPattern(pitch, p_geom.pupdiam * 1.1, pupAngle)
390 if p_dm.margin_out
is not None:
391 print(f
'p_dm.margin_out={p_dm.margin_out} is being '
392 'used for pupil-based actuator filtering')
393 pup_side = p_geom._ipupil.shape[0]
394 cub_off = dm_util.filterActuWithPupil(cub + pup_side // 2 - 0.5,
396 p_dm.margin_out * p_dm.get_pitch())
397 cub = cub_off - pup_side // 2 + 0.5
398 p_dm.set_ntotact(cub.shape[1])
399 elif p_dm.type_pattern == scons.PatternType.SQUARE:
400 print(
"Pattern type : square")
401 cub = dm_util.createSquarePattern(pitch, nxact + 4)
403 raise ValueError(
"This pattern does not exist for pzt dm")
406 inbigcirc = np.arange(cub.shape[1])
410 inbigcirc = dm_util.select_actuators(cub[0, :], cub[1, :], p_dm.nact,
411 p_dm._pitch, cobs, p_dm.margin_in,
412 p_dm.margin_out, p_dm._ntotact)
413 p_dm._ntotact = inbigcirc.size
421 cubval = cub[:, inbigcirc]
422 ntotact = cubval.shape[1]
426 i1t = (cubval[0, :] - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
427 j1t = (cubval[1, :] - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
436 influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
439 print(
"Computing Influence Function type : ", p_dm.influ_type)
441 for i
in tqdm(range(ntotact)):
444 x = np.tile(np.arange(i1, i1 + smallsize, dtype=np.float32),
452 y = np.tile(np.arange(j1, j1 + smallsize, dtype=np.float32),
458 if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
459 influ[:, :, i] = influ_util.makeRadialSchwartz(pitch, coupling, x=x, y=y)
460 elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
461 influ[:, :, i] = influ_util.makeSquareSchwartz(pitch, coupling, x=x, y=y)
462 elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
463 influ[:, :, i] = influ_util.makeBlacknutt(pitch, coupling, x=x, y=y)
464 elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
465 influ[:, :, i] = influ_util.makeGaussian(pitch, coupling, x=x, y=y)
466 elif (p_dm.influ_type == scons.InfluType.BESSEL):
467 influ[:, :, i] = influ_util.makeBessel(pitch, coupling, x=x, y=y,
468 patternType=p_dm.type_pattern)
469 elif (p_dm.influ_type == scons.InfluType.DEFAULT):
470 influ[:, :, i] = influ_util.makeRigaut(pitch, coupling, x=x, y=y)
472 print(
"ERROR influtype not recognized (defaut or gaussian or bessel)")
474 if (p_dm._puppixoffset
is not None):
475 xpos += p_dm._puppixoffset[0]
476 ypos += p_dm._puppixoffset[1]
478 if p_dm.segmented_mirror:
481 s = (p_geom._ipupil.shape[0] - p_geom._mpupil.shape[0]) // 2
483 from skimage.morphology
import label
485 for i
in tqdm(range(ntotact)):
487 i1, j1 = i1t[i] + s - smallsize // 2, j1t[i] + s - smallsize // 2
488 pupilSnapshot = p_geom._ipupil[i1:i1 + smallsize, j1:j1 + smallsize]
489 if np.all(pupilSnapshot):
491 labels, num = label(pupilSnapshot, background=0, return_num=
True)
495 maxPerArea = np.array([
496 (influ[:, :, i] * (labels == k).astype(np.float32)).max()
497 for k
in range(1, num + 1)
499 influ[:, :, i] *= (labels == np.argmax(maxPerArea) + 1).astype(np.float32)
500 print(f
'{k} cross-spider influence functions trimmed.')
502 influ = influ * float(p_dm.unitpervolt / np.max(influ))
508 dim = max(p_geom._mpupil.shape[0], p_dm._n2 - p_dm._n1 + 1)
509 off = (dim - p_dm._influsize) // 2
512 def init_custom_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, diam: float):
513 """Read Fits for influence pzt fonction and form
516 p_dm: (Param_dm) : dm settings
518 p_geom: (Param_geom) : geom settings
520 diam: (float) : tel diameter
522 Conversion. There are sereval coordinate systems.
523 Some are coming from the input fits file, others from compass.
524 Those systems differ by scales and offsets.
525 Coord systems from the input fits file:
526 - they all have the same scale: the coordinates are expressed in
528 - one system is the single common frame where fits data are described [i]
529 - one is local to the minimap [l]
530 Coord systems from compass:
531 - they all have the same scale (different from fits one) expressed in
533 - one system is attached to ipupil (i=image, largest support) [i]
534 - one system is attached to mpupil (m=medium, medium support) [m]
535 - one system is local to minimap [l)]
537 Variables will be named using
538 f, c: define either fits or compass
541 from astropy.io
import fits
as pfits
544 hdul = pfits.open(p_dm.file_influ_hdf5)
545 print(
"Read influence function from fits file : ", p_dm.file_influ_hdf5)
547 xC = hdul[0].header[
'XCENTER']
548 yC = hdul[0].header[
'YCENTER']
549 pitchPix = hdul[0].header[
'PITCHPX']
550 pitchMeters = hdul[0].header[
'PITCHM']
551 hi_i1, hi_j1 = hdul[1].data
553 xpos, ypos = hdul[3].data
558 pitchPixCompass = pitchMeters / p_geom._pixsize
559 scaleToCompass = pitchPixCompass / pitchPix
565 offsetXToCompass = p_geom.cent - xC * scaleToCompass
566 offsetYToCompass = p_geom.cent - yC * scaleToCompass
569 iSize, jSize, ntotact = influ.shape
573 ess1 = np.ceil((hi_i1+iSize) * scaleToCompass + offsetXToCompass) \
574 - np.floor(hi_i1 * scaleToCompass + offsetXToCompass)
577 ess2 = np.ceil((hi_j1+jSize) * scaleToCompass + offsetYToCompass) \
578 - np.floor(hi_j1 * scaleToCompass + offsetYToCompass)
580 smallsize = np.maximum(ess1, ess2).astype(int)
583 p_dm._ntotact = ntotact
584 p_dm._influsize = np.int(smallsize)
585 p_dm._i1 = np.zeros(ntotact, dtype=np.int32)
586 p_dm._j1 = np.zeros(ntotact, dtype=np.int32)
587 p_dm._xpos = np.zeros(ntotact, dtype=np.float32)
588 p_dm._ypos = np.zeros(ntotact, dtype=np.float32)
589 p_dm._influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
592 for i
in range(ntotact):
594 hi_x = np.arange(iSize) + hi_i1[i]
595 hi_y = np.arange(jSize) + hi_j1[i]
598 ci_x = hi_x * scaleToCompass + offsetXToCompass
599 ci_y = hi_y * scaleToCompass + offsetYToCompass
603 i] * scaleToCompass + offsetXToCompass
604 ci_j1 = hi_j1[i] * scaleToCompass + offsetYToCompass
605 ci_i1 = np.floor(ci_i1).astype(np.int32)
606 ci_j1 = np.floor(ci_j1).astype(np.int32)
607 ci_xpix = ci_i1 + np.arange(smallsize)
608 ci_ypix = ci_j1 + np.arange(smallsize)
613 f = interpolate.interp2d(ci_y, ci_x, influ[:, :, i], kind=
'cubic')
614 temp = f(ci_ypix, ci_xpix) * p_dm.unitpervolt
616 p_dm._influ[:, :, i] = temp
622 tmp = p_geom.ssize - (np.max(p_dm._i1 + smallsize))
623 margin_i = np.minimum(tmp, np.min(p_dm._i1))
624 tmp = p_geom.ssize - (np.max(p_dm._j1 + smallsize))
625 margin_j = np.minimum(tmp, np.min(p_dm._j1))
627 p_dm._xpos = xpos * scaleToCompass + offsetXToCompass
628 p_dm._ypos = ypos * scaleToCompass + offsetYToCompass
630 p_dm._n1 = int(np.minimum(margin_i, margin_j))
631 p_dm._n2 = p_geom.ssize - 1 - p_dm._n1
639 def make_tiptilt_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
641 """Compute the influence functions for a tip-tilt DM
644 p_dm: (Param_dm) : dm settings
646 patchDiam: (int) : patchDiam for dm size
648 p_geom: (Param_geom) : geom settings
650 diam: (float) : telescope diameter
652 influ: (np.ndarray(dims=3,dtype=np.float64)) : cube of the IF
655 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
659 influ = dm_util.make_zernike(nzer + 1, dim, patchDiam, p_geom.cent - p_dm._n1 + 1,
660 p_geom.cent - p_dm._n1 + 1, 1)[:, :, 1:]
663 current = influ[dim // 2 - 1, dim // 2 - 1, 0] - \
664 influ[dim // 2 - 2, dim // 2 - 2, 0]
665 fact = p_dm.unitpervolt * diam / p_geom.pupdiam * 4.848 / current
668 p_dm._ntotact = influ.shape[2]
669 p_dm._influsize = influ.shape[0]
675 def make_kl_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
676 cobs: float) ->
None:
677 """Compute the influence function for a Karhunen-Loeve DM
680 p_dm: (Param_dm) : dm settings
682 patchDiam: (int) : patchDiam for dm size
684 p_geom: (Param_geom) : geom settings
686 cobs: (float) : telescope cobs
689 dim = p_geom._mpupil.shape[0]
691 print(
"KL type: ", p_dm.type_kl)
694 nr = np.long(5.0 * np.sqrt(52))
695 npp = np.long(10.0 * nr)
697 nr = np.long(5.0 * np.sqrt(p_dm.nkl))
698 npp = np.long(10.0 * nr)
700 radp = kl_util.make_radii(cobs, nr)
702 kers = kl_util.make_kernels(cobs, nr, radp, p_dm.type_kl, p_dm.outscl)
704 evals, nord, npo, ordd, rabas = kl_util.gkl_fcom(kers, cobs, p_dm.nkl)
706 azbas = kl_util.make_azimuth(nord, npp)
708 ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap = kl_util.set_pctr(
709 patchDiam, nr, npp, p_dm.nkl, cobs, nord)
711 p_dm._ntotact = p_dm.nkl
720 p_dm._i1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
721 (dim - patchDiam) // 2
722 p_dm._j1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
723 (dim - patchDiam) // 2
724 p_dm._ntotact = p_dm.nkl
728 def comp_dmgeom(p_dm: conf.Param_dm, p_geom: conf.Param_geom):
729 """Compute the geometry of a DM : positions of actuators and influence functions
732 dm: (Param_dm) : dm settings
734 geom: (Param_geom) : geom settings
736 smallsize = p_dm._influsize
738 dm_dim = int(p_dm._n2 - p_dm._n1 + 1)
739 mpup_dim = p_geom._mpupil.shape[0]
741 if (dm_dim < mpup_dim):
742 offs = (mpup_dim - dm_dim) // 2
747 indgen = np.tile(np.arange(smallsize, dtype=np.int32), (smallsize, 1))
749 tmpx = np.tile(indgen, (nact, 1, 1))
750 tmpy = np.tile(indgen.T, (nact, 1, 1))
752 tmpx += offs + p_dm._i1[:,
None,
None]
753 tmpy += offs + p_dm._j1[:,
None,
None]
755 tmp = tmpx + mpup_dim * tmpy
758 tmp[tmpx < 0] = mpup_dim * mpup_dim + 10
759 tmp[tmpy < 0] = mpup_dim * mpup_dim + 10
760 tmp[tmpx > dm_dim - 1] = mpup_dim * mpup_dim + 10
761 tmp[tmpy > dm_dim - 1] = mpup_dim * mpup_dim + 10
762 itmps = np.argsort(tmp.flatten()).astype(np.int32)
763 tmps = tmp.flatten()[itmps].astype(np.int32)
764 itmps = itmps[np.where(itmps > -1)]
766 istart = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
767 npts = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
769 tmps_unique, cpt = np.unique(tmps, return_counts=
True)
770 if (tmps_unique > npts.size - 1).any():
771 tmps_unique = tmps_unique[:-1]
774 for i
in range(tmps_unique.size):
775 npts[tmps_unique[i]] = cpt[i]
776 istart[1:] = np.cumsum(npts[:-1])
778 p_dm._influpos = itmps[:np.sum(npts)].astype(np.int32)
782 p_dm._ninflu = npts.astype(np.int32)
783 p_dm._influstart = istart.astype(np.int32)
795 def correct_dm(context, dms: Dms, p_dms: list, p_controller: conf.Param_controller,
796 p_geom: conf.Param_geom, imat: np.ndarray =
None, dataBase: dict = {},
797 use_DB: bool =
False):
798 """Correct the geometry of the DMs using the imat (filter unseen actuators)
801 context: (carmaWrap_context): context
802 dms: (Dms) : Dms object
803 p_dms: (list of Param_dm) : dms settings
804 p_controller: (Param_controller) : controller settings
805 p_geom: (Param_geom) : geom settings
806 imat: (np.ndarray) : interaction matrix
807 dataBase: (dict): dictionary containing paths to files to load
808 use_DB: (bool): dataBase use flag
810 print(
"Filtering unseen actuators... ")
812 resp = np.sqrt(np.sum(imat**2, axis=0))
814 ndm = p_controller.ndm.size
817 for nmc
in range(ndm):
818 nm = p_controller.ndm[nmc]
819 nactu_nm = p_dms[nm]._ntotact
821 if (p_dms[nm].type == scons.DmType.PZT):
823 influpos, ninflu, influstart, i1, j1, ok = h5u.load_dm_geom_from_dataBase(
825 p_dms[nm].set_ntotact(ok.shape[0])
826 p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
827 p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
828 p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
829 p_dms[nm]._influpos = influpos
830 p_dms[nm]._ninflu = ninflu
831 p_dms[nm]._influstart = influstart
835 tmp = resp[inds:inds + p_dms[nm]._ntotact]
836 ok = np.where(tmp > p_dms[nm].thresh * np.max(tmp))[0]
837 nok = np.where(tmp <= p_dms[nm].thresh * np.max(tmp))[0]
839 p_dms[nm].set_ntotact(ok.shape[0])
840 p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
841 p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
842 p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
843 p_dms[nm].set_i1(p_dms[nm]._i1[ok])
844 p_dms[nm].set_j1(p_dms[nm]._j1[ok])
848 h5u.save_dm_geom_in_dataBase(nmc, p_dms[nm]._influpos,
850 p_dms[nm]._influstart, p_dms[nm]._i1,
853 dim = max(p_dms[nm]._n2 - p_dms[nm]._n1 + 1, p_geom._mpupil.shape[0])
854 ninflupos = p_dms[nm]._influpos.size
855 n_npts = p_dms[nm]._ninflu.size
857 dms.insert_dm(context, p_dms[nm].type, p_dms[nm].alt, dim,
858 p_dms[nm]._ntotact, p_dms[nm]._influsize, ninflupos, n_npts,
859 p_dms[nm].push4imat, 0, p_dms[nm].dx / p_geom._pixsize,
860 p_dms[nm].dy / p_geom._pixsize, p_dms[nm].theta, p_dms[nm].G,
861 context.active_device, nm)
862 dms.d_dms[nm].pzt_loadarrays(p_dms[nm]._influ, p_dms[nm]._influpos.astype(
863 np.int32), p_dms[nm]._ninflu, p_dms[nm]._influstart, p_dms[nm]._i1,
872 makePetalDm(p_dm, p_geom, pupAngleDegree)
874 The function builds a DM, segmented in petals according to the pupil
875 shape. The petals will be adapted to the EELT case only.
877 <p_geom> : compass object p_geom. The function requires the object p_geom
878 in order to know what is the pupil mask, and what is the mpupil.
879 <p_dm> : compass petal dm object p_dm to be created. The function will
880 transform/modify in place the attributes of the object p_dm.
884 p_dm._n1 = p_geom._n1
885 p_dm._n2 = p_geom._n2
887 p_dm._influsize = smallsize
888 p_dm.set_ntotact(nbSeg)
891 p_dm._xpos = i1 + smallsize / 2 + p_dm._n1
892 p_dm._ypos = j1 + smallsize / 2 + p_dm._n1
901 <pupImage> : image of the pupil
903 La fonction renvoie des fn d'influence en forme de petale d'apres
904 une image de la pupille, qui est supposee etre segmentee.
907 influ, i1, j1, smallsize, nbSeg = make_petal_dm_core(pupImage, 0.0)
913 from scipy.ndimage.measurements
import label
914 from scipy.ndimage.morphology
import binary_opening
915 s = np.ones((2, 2), dtype=np.bool)
916 segments, nbSeg = label(binary_opening(pupImage, s))
925 for i
in range(nbSeg):
926 petal = segments == (i + 1)
927 profil = np.sum(petal, axis=1) != 0
928 extent = np.sum(profil).astype(np.int32)
929 i1t.append(np.min(np.where(profil)[0]))
930 i2t.append(np.max(np.where(profil)[0]))
931 if extent > smallsize:
934 profil = np.sum(petal, axis=0) != 0
935 extent = np.sum(profil).astype(np.int32)
936 j1t.append(np.min(np.where(profil)[0]))
937 j2t.append(np.max(np.where(profil)[0]))
938 if extent > smallsize:
945 influ = np.zeros((smallsize, smallsize, nbSeg), dtype=np.float32)
947 npt = pupImage.shape[0]
948 i0 = j0 = npt / 2 - 0.5
949 petalMap =
build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
950 ii1 = np.zeros(nbSeg)
951 jj1 = np.zeros(nbSeg)
952 for i
in range(nbSeg):
953 ip = (smallsize - i2t[i] + i1t[i] - 1) // 2
954 jp = (smallsize - j2t[i] + j1t[i] - 1) // 2
955 i1 = np.maximum(i1t[i] - ip, 0)
956 j1 = np.maximum(j1t[i] - jp, 0)
957 if (j1 + smallsize) > npt:
959 if (i1 + smallsize) > npt:
962 k = petalMap[i1 + smallsize // 2, j1 + smallsize // 2]
963 petal = (petalMap == k)
964 influ[:, :, k] = petal[i1:i1 + smallsize, j1:j1 + smallsize]
968 return influ, ii1, jj1, int(smallsize), nbSeg
973 Makes an image npt x npt of <nbSeg> regularly spaced angular segments
975 Origin of angles is set by <pupAngleDegree>.
977 The segments are oriented as defined in document "Standard Coordinates
978 and Basic Conventions", ESO-193058.
979 This document states that the X axis lies in the middle of a petal, i.e.
980 that the axis Y is along the spider.
981 The separation angle between segments are [-30, 30, 90, 150, -150, -90].
982 For this reason, an <esoOffsetAngle> = -pi/6 is introduced in the code.
988 p = build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
991 rot = pupAngleDegree * np.pi / 180.0
994 esoOffsetAngle = -np.pi / 6
995 x = np.arange(npt) - i0
996 y = np.arange(npt) - j0
997 X, Y = np.meshgrid(x, y, indexing=
'ij')
998 theta = (np.arctan2(Y, X) - rot + 2 * np.pi - esoOffsetAngle) % (2 * np.pi)
1001 angleStep = 2 * np.pi / nbSeg
1002 startAngle = np.arange(nbSeg) * angleStep
1003 endAngle = np.roll(startAngle, -1)
1004 endAngle[-1] = 2 * np.pi
1005 petalMap = np.zeros((npt, npt), dtype=int)
1006 for i
in range(nbSeg):
1007 nn = np.where(np.logical_and(theta >= startAngle[i], theta < endAngle[i]))